!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!
!
#ifdef VAR_DEBUG
#define HSODEBUG .TRUE.
#else
#define HSODEBUG .FALSE.
#endif
C  /* Deck hsoctl */
      SUBROUTINE HSOCTL (WORD,GP,CMO,UDV,PV,XINDX,ANTSYM,WRK,LWRK)
C
C Olav Vahtras
C Jun 20, 1990
C
C Driver routine for the construction of spin-orbit property vector
C
C Input: Spin-orbit component (as read from LUPROP), WORD
C        MO coefficients,CMO
C        First order reduced density matrix, UDV
C        Second order reduced density matrix, PV
C        Index vector, XINDX
C
C Output: Spin-orbit property vector returned in first elements of WRK
C
#include "implicit.h"
      DIMENSION  GP(*),CMO(*),UDV(*),PV(*),XINDX(*),WRK(LWRK)
      CHARACTER*8 WORD
C
      PARAMETER (D0=0.D0, D1=1.D0, DM1=-1.D0, D2=2.D0)
      LOGICAL    SO2TRA, FILE_EXISTS, FOPEN, DO1,DO2, LORB, LCON, TRPDEN
      CHARACTER*8 LABEL
      LOGICAL    FNDLAB
      EXTERNAL   FNDLAB
C
C Used from common blocks:
C  MULD2H,ISMO,ISW,IOBTYP,IDBTYP,NISHT,NORBT,NASHT,NNASHX
C  INFORB: MULD2H,NISHT,NASHT,NORBT,NNASHX
C  INFIND: ISW,IOBTYP
C  TRHSO : ILXYZ, KSYMSO, OLDTRA
C  INFHSO: IPRHSO, TESTZY, DOSO1, DOSO2
C  INFHYP: HYPCAL
C  INFSMO: SOMOM
C  INFPP : EXMOM
C  INFRSP: IPRRSP,SOPPA,???
C
C-- common blocks:
#include "maxash.h"
#include "maxorb.h"
#include "aovec.h"
#include "priunit.h"
#include "dummy.h"
#include "inforb.h"
#include "inftap.h"
#include "infind.h"
#include "infvar.h"
#include "infrsp.h"
#include "wrkrsp.h"
#include "infpri.h"
#include "trhso.h"
#include "rspprp.h"
#include "infhso.h"
#include "infhyp.h"
#include "infsmo.h"
#include "infpp.h"
#include "iratdef.h"
#include "codata.h"
#include "inflr.h"
#include "qrinf.h"
#include "cbihr2.h"
#include "infesr.h"
#include "dftcom.h"
#include "infrank.h"
C
      CALL QENTER('HSOCTL')
C
      IF (SOPPA) CALL QUIT('HSOCTL: SOPPA not implemented yet!')
      IPRHSO = MAX(IPRHSO,IPRRSP,IPRPP,IPRLR,IPRESR)
      KSYMSO=KSYMOP
      CALL DZERO(GP,KZYVAR)
      HSOFAC=ALPHAC**2/4
      IF (WORD(2:2).EQ.'1') THEN
         OPRANK(INDPRP(WORD)) =1
         CALL QRGP(WORD,GP,CMO,XINDX,ANTSYM,WRK,LWRK)
         CALL DSCAL(KZYVAR,HSOFAC,GP,1)
         GO TO 9999
      END IF
      ANTSYM = 1.0D0
      IF (WORD(2:2).EQ.'2') THEN
         DO1 = .FALSE.
         DO2 = .TRUE.
      END IF
      IF (WORD(2:2).EQ.' ') THEN
         DO1 = DOSO1
         DO2 = DOSO2
      END IF
C
      IF (IPRHSO.GT.2) THEN
         CALL TIMER('START ',HSOSTA,HSOTIM)
         CALL HEADER('Output from HSOCTL',-1)
         WRITE(LUPRI,'(/2A,3X,A)')
     *   ' Spin-orbit property vector calculation',
     *   ' component = ',WORD
         WRITE(LUPRI,'(/A,I5)')' Print level in HSOCTL: ',IPRHSO
         IF (TESTZY) WRITE(LUPRI,'(/A)')
     *' Z and Y parts of configurational property vector explicitly'
         IF (.NOT.DO1) WRITE(LUPRI,'(/A)')
     *' Skip one-electron spin-orbit contributions'
         IF (.NOT.DO2) WRITE(LUPRI,'(/A)')
     *' Skip two-electron spin-orbit contributions'
      END IF
      LORB = KZWOPT.GT.0
      IF (KSYMOP.EQ.1) THEN
         LCON = KZCONF.GT.1
      ELSE
         LCON = KZCONF.GT.0
      END IF
C
C Check if gradient is on file
C
      INQUIRE(FILE='HSOGRAD',EXIST=FILE_EXISTS)
      LUHSO = -1
      CALL GPOPEN(LUHSO,'HSOGRAD','UNKNOWN',' ','UNFORMATTED',IDUMMY,
     &            .FALSE.)
      REWIND LUHSO
      IF (FILE_EXISTS) THEN
C     ... aug07-hjaaj: gfortran doesn't like BACKSPACE(LUHSO) in FNDLAB on empty file!
      IF (FNDLAB(WORD,LUHSO)) THEN
         CALL READT(LUHSO,KZYVAR,GP)
         CALL GPCLOSE(LUHSO,'KEEP')
         GO TO 9999
      END IF
      END IF
C
C ALLOCATE WORK SPACE
C
      KWRK1 = 1
      LWRK1 = LWRK
      CALL MEMGET2('REAL','FC', KFC    , NORBT*NORBT,WRK,KWRK1,LWRK1)
      IF (LORB) THEN
         CALL MEMGET2('REAL','FV', KFV    , NORBT*NORBT,WRK,KWRK1,LWRK1)
         CALL MEMGET2('REAL','QA', KQA    , NORBT*NASHT,WRK,KWRK1,LWRK1)
         CALL MEMGET2('REAL','QB', KQB    , NORBT*NASHT,WRK,KWRK1,LWRK1)
      END IF
      CALL MEMGET2('REAL','H1', KH1 , NORBT*NORBT,WRK,KWRK1,LWRK1)
      CALL MEMGET2('REAL','H2', KH2 , NORBT*NORBT,WRK,KWRK1,LWRK1)
      IF (LCON) THEN
         CALL MEMGET2('REAL','H2AC',KH2AC,NNASHX*NNASHX,WRK,KWRK1,LWRK1)
      ELSE
         CALL MEMGET2('REAL','H2AC',KH2AC,0            ,WRK,KWRK1,LWRK1)
      END IF
C
      IF (IPRHSO.GT.20 .AND. LORB .AND. NASHT.GT.0) THEN
            WRITE(LUPRI,'(/A)')
     *   ' First order density matrix:'
            CALL OUTPUT(UDV,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI)
            WRITE(LUPRI,'(/A)')
     *   ' Second order density matrix:'
            CALL PRIAC2(PV,NASHT,LUPRI)
      END IF
C
C
C Read AO one-electron property integrals and transform to MO basis.
C
      IF (DO1) THEN
         LABEL = WORD
         LABEL(2:2) = '1'
         KSYMP = -1
         CALL PRPGET (LABEL,CMO,WRK(KH1),KSYMP,ANTSYM,WRK(KWRK1),LWRK1,
     &                IPRHSO)
         IF (KSYMP.NE.KSYMOP) THEN
            WRITE (LUPRI,'(/A/2A/A,2I5/A,F10.2)')
     &           'FATAL ERROR: KSYMOP .ne. KSYMP from PRPGET',
     &           '   Property label  : ',WORD,
     &           '   KSYMOP and KSYMP:',KSYMOP,KSYMP,
     &           '   ANTSYM          :',ANTSYM
            CALL QUIT('KSYMOP .ne. KSYMP from PRPGET')
         END IF
      END IF
C
C Print atomic and molecular property integrals, if desired
C
      IF (IPRHSO.GT.20 .AND. DO1) THEN
         WRITE (LUPRI,'(/2A)')' Atomic property integrals:', LABEL
         CALL OUTPUT(WRK(KWRK1),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
         WRITE (LUPRI,'(/2A)')' Molecular property integrals:', LABEL
         CALL OUTPUT(WRK(KH1),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
      END IF
C     IF (.NOT.DO2) GOTO 95
      IF (.NOT.TDHF) CALL TRANSFORM_HSO(WORD,CMO,WRK,KWRK1,LWRK1)
C
      IF (TRPLET) THEN
         TRPDEN=.FALSE.
         CALL HSOFCK(WORD,WRK(KFC),WRK(KFV),WRK(KQA),WRK(KQB),
     &               WRK(KH2AC),UDV,PV,CMO,LORB,LCON,WRK,KWRK1,LWRK1)
      ELSE
         TRPDEN=.TRUE.
         CALL MEMGET2('REAL','D',KD,N2ASHX,WRK,KWRK1,LWRK1)
         CALL MEMGET2('REAL','P',KP,2*N2ASHX*N2ASHX,WRK,KWRK1,LWRK1)
         CALL DZERO(WRK(KD),N2ASHX)
         CALL DZERO(WRK(KP),2*N2ASHX*N2ASHX)
         IF (TDHF) THEN
            IF (NASHT.GE.1) THEN
               CALL DUNIT(WRK(KD),NASHT)
            END IF
         ELSE
C
C Need new density matrices
C
            CALL MEMGET2('REAL','CREF',KCREF,MZCONF(1),WRK,KWRK1,LWRK1)
            CALL GETREF(WRK(KCREF),MZCONF(1))
            CALL RSPDM(IREFSY,IREFSY,MZCONF(1),MZCONF(1),
     &         WRK(KCREF),WRK(KCREF),
     &         WRK(KD),WRK(KP),1,0,.FALSE.,.FALSE.,XINDX,
     &         WRK,KWRK1,LWRK1)
            KP1=KP
            KP2=KP+N2ASHX*N2ASHX
            CALL MTRSP(N2ASHX,N2ASHX,WRK(KP1),N2ASHX,WRK(KP2),N2ASHX)
            CALL MEMREL('HSOCTL<-RSPDM',WRK,KD,KCREF,KWRK1,LWRK1)
         END IF
         CALL HSOFCK(WORD,WRK(KFC),WRK(KFV),WRK(KQA),WRK(KQB),
     &   WRK(KH2AC),WRK(KD),WRK(KP),CMO,LORB,LCON,WRK,KWRK1,LWRK1)
         CALL MEMREL('HSOCTL<-HSOFCK',WRK,KD,KP,KWRK1,LWRK1)
      END IF
 95   CONTINUE
C
C
C Add the one-electron and two electron parts of the inactive Fock matrix
C
      IF (DO1) THEN
         CALL DAXPY(N2ORBX,D1,WRK(KH1),1,WRK(KFC),1)
         IF (TDHF.AND.NASHT.GT.0) THEN
            CALL DAXPY(N2ORBX,DM1,WRK(KH1),1,WRK(KFV),1)
         END IF
      END IF
C
C
C
C Construct orbital property vector
C
      IF (LORB) THEN
         IF (IPRHSO.GT.2) CALL TIMER('START ',ORBSTA,ORBTIM)
         IF (TRPLET) THEN
            CALL HSOORB(.TRUE.,1,DUMMY1,
     *         WRK(KFC),WRK(KFV),WRK(KQA),WRK(KQB),
     *         UDV    ,GP,TRPDEN)
         ELSE
C
C Quick fix for high spin matrices which do not fit neatly into HSOORB
C
            IF (TDHF.AND.NASHT.GT.0) TRPDEN=.FALSE.
C
            CALL HSOORB(.TRUE.,1,DUMMY1,
     *         WRK(KFC),WRK(KFV),WRK(KQA),WRK(KQB),
     *         WRK(KD),GP,TRPDEN)
         END IF
C        CALL HSOORB(ONEIND,NSIM,FC,FCX,FVX,QAX,QBX,UDV,
C    *                  EVECS,TRPDEN)
C
C
         IF (IPRHSO.GT.2) CALL TIMER('HSOORB',ORBSTA,ORBTIM)
      END IF
C
C  Compensate for the sign in HSOORB
C
      CALL DSCAL(KZWOPT,DM1,GP(1+KZCONF),1)
      CALL DSCAL(KZWOPT,DM1,GP(1+KZCONF+KZVAR),1)
C
C Construct configuration property vector
C
      IF (LCON) THEN
         IF (IPRHSO.GT.2) CALL TIMER('START ',ORBSTA,ORBTIM)
         CALL HSOSIG(WRK(KFC),WRK(KH2AC),
     *      GP,XINDX,WRK(KWRK1),LWRK1)
C
C        CALL HSOSIG(FC,H2AC, GP,XINDX,WRK,LWRK)
C
C
         IF (IPRHSO.GT.2) CALL TIMER('HSOSIG',ORBSTA,ORBTIM)
      END IF
      IF (IPRHSO.GT.10) THEN
         WRITE(LUPRI,'(//A//A)')
     *      ' Configuration property vector: ',
     *      '               Z part         Y part'
         CALL OUTPUT(GP,1,KZCONF,1,2,KZVAR,2,1,LUPRI)
         WRITE(LUPRI,'(//A//A)')
     *      ' Orbital property vector: ',
     *      '               Z part         Y part'
         CALL OUTPUT(GP(1+KZCONF),1,KZWOPT,1,2,KZVAR,2,1,LUPRI)
      ELSE IF (IPRHSO.GT.2) THEN
         IF (LORB) CALL RSPPRO (GP(1+KZCONF),KZVAR,UDV,LUPRI)
         IF (LCON) CALL RSPPRC (GP,KZCONF,KZVAR,LUPRI)
      END IF
C
      IF (IPRHSO.GT.2) CALL TIMER('HSOCTL',HSOSTA,HSOTIM)
      IF (X2GRAD) THEN
         CALL MEMREL('X2GRAD',WRK,KFC,KFC,KWRK1,LWRK1)
         CALL MEMGET2('INTE','MJWOP',KMJWOP,16*MAXWOP,WRK,KWRK1,LWRK1)
         CALL HEADER('X2GRAD TEST FOR SPIN-ORBIT GRADIENT ELEMENTS',3)
         CALL MEMGET2('REAL','GP2',KGP2,KZYVAR,WRK,KWRK1,LWRK1)
         CALL DZERO(WRK(KGP2),KZYVAR)
         CALL SETZY(WRK(KMJWOP))
         CALL HSOAL2 (WORD,WRK(KGP2),CMO,UDV,PV,XINDX,WRK(KMJWOP),
     &                WRK(KWRK1),LWRK1)
         WRITE(LUPRI,'(/A)')'               GP1           GP2'
         WRITE(LUPRI,'(/A)')'      ---------------------------------'
         DMAXGP = D0
         DO 77, K = 0,KZYVAR-1
            WRITE(LUPRI,'(10X,2F14.8)') GP(1+K), WRK(KGP2+K)
            DMAXGP = MAX(DMAXGP,ABS(GP(1+K)-WRK(KGP2+K)))
77       CONTINUE
         WRITE(LUPRI,'(/A)')'      ---------------------------------'
         WRITE(LUPRI,'(/A,E20.8)')
     *      'LARGEST DIFFERENCE OF SPIN-ORBIT GRADIENT ELEMENTS' ,
     *      DMAXGP
         CALL MEMREL('X2GRAD',WRK,KGP2,KGP2,KWRK1,LWRK1)
      END IF
      CALL DSCAL(KZYVAR,HSOFAC,GP,1)
C
C Save on file
C
      REWIND LUHSO
      CALL NEWLAB(WORD,LUHSO,LUERR)
      CALL WRITT(LUHSO,KZYVAR,GP)
      CALL MEMREL('HSOCTL',WRK,1,1,KWRK1,LWRK1)
      CALL GPCLOSE(LUHSO,'KEEP')
C
 9999 CALL QEXIT('HSOCTL')
      RETURN
      END
C  /* Deck hsoinp */
      SUBROUTINE HSOINP(WORD)
C
#include "implicit.h"
C
#include "priunit.h"
#include "infrsp.h"
#include "wrkrsp.h"
#include "infs0.h"
#include "infpri.h"
#include "infdim.h"
#include "trhso.h"
#include "infhso.h"
C
      LOGICAL NEWDEF
      PARAMETER ( NTABLE =  8 )
      CHARACTER PROMPT*1, WORD*7, TABLE(NTABLE)*7, WORD1*7
C
      DATA TABLE /'.TESTZY', '.SO1ONL', '.SO2ONL', '.PRINT', 
     *            '.OLDTRA', '.X2MAT ', '.A2MAT ', '.X2GRAD'/
C
      NEWDEF = (WORD .EQ. '*SPIN-O')
      ICHANG = 0
      IF (NEWDEF) THEN
         WORD1 = WORD
  100    CONTINUE
            READ (LUCMD, '(A7)') WORD
            CALL UPCASE(WORD)
            PROMPT = WORD(1:1)
            IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') GO TO 100
            IF (PROMPT .EQ. '.') THEN
               ICHANG = ICHANG + 1
               DO I=1, NTABLE
                  IF (TABLE(I) .EQ. WORD) THEN
                     GO TO (1,2,3,4,5,6,7,8), I
                  END IF
               END DO
               IF (WORD .EQ. '.OPTION') THEN
                 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI)
                 GO TO 100
               END IF
               WRITE (LUPRI,'(/,3A,/)') ' KEYWORD "',WORD,
     *            '" NOT RECOGNIZED IN HSOINP.'
               CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI)
               CALL QUIT(' ILLEGAL KEYWORD IN HSOINP ')
 1             CONTINUE
                  TESTZY = .TRUE.
               GO TO 100
 2             CONTINUE
                  DOSO2 = .FALSE.
               GO TO 100
 3             CONTINUE
                  DOSO1 = .FALSE.
               GO TO 100
 4             CONTINUE
                  READ(LUCMD,*)IPRHSO
               GO TO 100
 5             CONTINUE
                  OLDTRA = .TRUE.
               GO TO 100
 6             CONTINUE
                  X2MAT = .TRUE.
               GO TO 100
 7             CONTINUE
                  A2MAT = .TRUE.
               GO TO 100
 8             CONTINUE
                  X2GRAD = .TRUE.
               GO TO 100
            ELSE IF (PROMPT.EQ.'*') THEN
               GO TO 300
            ELSE
               WRITE (LUPRI,'(/,3A,/)') ' PROMPT "',WORD,
     *            '" NOT RECOGNIZED IN RSPINP.'
               CALL QUIT(' ILLEGAL PROMPT IN HSOINP ')
            END IF
         GO TO 100
      END IF
  300 CONTINUE
      IF (ICHANG .GT. 0) THEN
         CALL HEADER('CHANGES OF DEFAULTS FOR HSOINP:',0)
         IF (TESTZY) WRITE(LUPRI,'(/A/A,L1)')
     *      ' Both parts of configuration property vector explicitely',
     *      ' TESTZY = ',TESTZY
         IF (IPRHSO.NE.2) WRITE(LUPRI,'(/A/A,I5)')
     *      ' Print level in spin-orbit property vector calculation',
     *      ' IPRHSO = ',IPRHSO
         IF (OLDTRA)  WRITE(LUPRI,'(/A/A,L1)')
     *      ' Use existing two-electron spin-orbit integral file',
     *      ' OLDTRA = ',OLDTRA
         IF (.NOT.DOSO1) WRITE(LUPRI,'(/A)')
     *      'Skip one-particle part in HSOCTL'
         IF (.NOT.DOSO2) WRITE(LUPRI,'(/A)')
     *      'Skip two-particle part in HSOCTL'
         IF (X2MAT) WRITE(LUPRI,'(/A)')
     *      'X2MAT: Calculate full X2 matrix (quadratic response)'
         IF (A2MAT) WRITE(LUPRI,'(/A)')
     *      'A2MAT: not implemented'
      END IF
C
C *** END OF HSOINP
C
      RETURN
      END
C  /* Deck fsomu */
      SUBROUTINE FSOMU(ICI1,IDI1,H2,FCSO,FVSO,UDV,LORB,WRK,LWRK)
C
C Olav Vahtras
C Apr 11, 1990
C
C CALCULATE ALL CONTRIBUTIONS TO INACTIVE AND ACTIVE SPIN-ORBIT
C FOCK MATRICES FROM MULLIKEN DISTRIBUTIONS
C
C  FCSO(P,Q) = SUM (K) 2*(KK|P^Q) - SUM (K) 3*(PK|K^Q)
C                                 - SUM (K) 3*(KQ|P^K)
C
C  FVSO(P,Q) = SUM (X,Y) (XY|P^Q) D(XY) - SUM (K) 3/2*(PX|Y^Q) D(XY)
C                                     - SUM (K) 3/2*(XQ|P^Y) D(XY)
C
C
#include "implicit.h"
C
C Used from common blocks:
C
C   INFORB : NORBT,NISHT,NASHT,MULD2H
C   INFIND : ISMO,IOBTYP,ICH,ISMO,IASH,IORB,NISH,NASH,NOCC,NORB
C   INFHSO :
C   WRKRSP : KSYMOP
C
#include "maxash.h"
#include "maxorb.h"
#include "priunit.h"
#include "inforb.h"
#include "infind.h"
#include "wrkrsp.h"
#include "infrsp.h"
#include "infpri.h"
#include "orbtypdef.h"
#include "trhso.h"
#include "infhso.h"
C
      DIMENSION H2(NORBT,NORBT)
      DIMENSION FCSO(NORBT,NORBT),FVSO(NORBT,NORBT)
      DIMENSION UDV(NASHT,NASHT),WRK(*)
      LOGICAL LORB
C
      PARAMETER (D1=1.0D0, D1P5=1.5D0, D2=2.0D0)
C
      CALL QENTER('FSOMU')
C
C
C  Local print level
C
      IPRFSO = 20
C
C     Order (C,D) index such that C .ge. D
C     in inactive-active-secondary order (using ISW)
C
      IF (ISW(ICI1) .GE. ISW(IDI1)) THEN
         ICI = ICI1
         IDI = IDI1
         DISFAC=D1
      ELSE
         ICI = IDI1
         IDI = ICI1
         DISFAC=-D1
      END IF
      IF (TRPLET) THEN
         COUFAC=D1
      ELSE
         COUFAC=D2
      END IF
C
C     Find distribution type ITYPCD
C
      ITYPC  = IOBTYP(ICI)
      ITYPD  = IOBTYP(IDI)
      ITYPCD = IDBTYP(ITYPC,ITYPD)
C
C     We only need secondary-occupied distributions
C
      IF (ITYPCD .EQ. JTSESE) GO TO 9999
C
      IF (ITYPC .EQ. JTACT) NCIW = ICH(ICI)
      IF (ITYPD .EQ. JTACT) NDIW = ICH(IDI)
C
      ICSYM = ISMO(ICI)
      IDSYM = ISMO(IDI)
      ICDSYM = MULD2H(ICSYM,IDSYM)
      KCDSYM = MULD2H(KSYMOP,ICDSYM)
C
      IF (IPRHSO.GT.IPRFSO) THEN
         WRITE(LUPRI,'(/A)')' ------ Output from FSOMU ------'
         WRITE(LUPRI,'(/A,2I5,5X,A,2X,A)')' Distribution CD',ICI1,IDI1,
     *                          COBTYP(ITYPC),COBTYP(ITYPD)
         WRITE(LUPRI,'(A,2I5)')' Reordered      ',ICI,IDI
         WRITE(LUPRI,'(A,2I5)')' Symmetry       ',ICSYM,IDSYM
      ENDIF
C
C
C Inactive Fock matrix
C
C Direct terms
C
C  FCSO(P,Q) = FCSO(P,Q) + SUM(K) 2*(KK|P^Q)
C
C here:
C        SUM(K) 2*(KK|C^D) -> FCSO(C,D)
C      - SUM(K) 2*(KK|C^D) -> FCSO(D,C)
C
         IF ( KSYMOP.EQ.ICDSYM ) THEN
            DO 10 I=1,NISHT
               IX=ISX(I)
               WRK(I)=H2(IX,IX)
10          CONTINUE
               IF (IPRHSO.GT.IPRFSO) THEN
                  WRITE(LUPRI,'(/A)')' Inactive direct terms'
                  WRITE(LUPRI,'(A)')' (KK|CD) diagonal'
                  CALL OUTPUT(WRK(1),1,NISHT,1,1,NISHT,1,1,LUPRI)
               END IF
               FAC=2*DISFAC*DSUM(NISHT,WRK(1),1)
               FCSO(ICI,IDI)=FCSO(ICI,IDI) + FAC
               FCSO(IDI,ICI)=FCSO(IDI,ICI) - FAC
         ENDIF
C
C Exchange terms: rearranged with D position inactive
C
C  FCSO(P,Q) = FCSO(P,Q) + SUM(K) 3*(PK|Q^K) - 3*(QK|P^K)
C
C
         IF (IPRHSO.GT.IPRFSO) THEN
            WRITE(LUPRI,'(/A)')' Loop over symmetries'
         END IF
         DO 200 ISYM = 1,NSYM
            IPSYM = ISYM
            IOFFP = IORB(IPSYM)
            NASHP = NASH(IPSYM)
            NOCCP = NOCC(IPSYM)
            NORBP = NORB(IPSYM)
            IASHP=IASH(IPSYM)
            IOFFPA=IOFFP+NISH(IPSYM)
            IF (IPRHSO.GT.IPRFSO) THEN
               WRITE(LUPRI,'(/A,I5)')' IPSYM',IPSYM
            ENDIF
            IF ( (NORBP.EQ.0) ) GO TO 200
C
            ICPSYM = MULD2H(ICSYM,IPSYM)
            IDPSYM = MULD2H(IDSYM,IPSYM)
C
C For the case D inactive
C
            IF ((ITYPCD.EQ.JTININ).OR.
     *          (ITYPCD.EQ.JTACIN).OR.
     *          (ITYPCD.EQ.JTSEIN)) THEN
C
C here:
C         + 3*(PD|C^D) -> FCSO(P,C)
C         - 3*(PD|C^D) -> FCSO(C,P)
C
C such that (PC) is at most secondary-active
C
               IF (IPRHSO.GT.IPRFSO) THEN
                  IF (ICPSYM.EQ.KSYMOP .AND.
     *                (ITYPCD.EQ.JTSEIN .AND. NOCCP.GT.0
     *                        .OR.
     *                      ITYPCD.NE.JTSEIN )
     *                        .OR.
     *                IDPSYM.EQ.KSYMOP .AND. ITYPCD.EQ.JTININ) THEN
                     WRITE(LUPRI,'(/A)')' Inactive exchange terms'
                  END IF
               END IF
               FAC = 3*DISFAC
               IF ( ICPSYM.EQ.KSYMOP ) THEN
                  IF (ITYPCD.EQ.JTSEIN) THEN
                     NDIMP = NOCCP
                  ELSE
                     NDIMP = NORBP
                  ENDIF
                  IF (NDIMP.GT.0) THEN
                     CALL DAXPY(NDIMP,FAC,H2(IOFFP+1,IDI),1,
     *                                 FCSO(IOFFP+1,ICI),1)
                     CALL DAXPY(NDIMP,-FAC,H2(IOFFP+1,IDI),1,
     *                                 FCSO(ICI,IOFFP+1),NORBT)
                     IF (IPRHSO.GT.IPRFSO) THEN
                       WRITE(LUPRI,'(A,I3,A,I3,I5)')' PC contribution ',
     *                    IOFFP+1,':',IOFFP+NDIMP,ICI
                       CALL OUTPUT(H2(IOFFP+1,IDI),1,NDIMP,1,1,
     *                                             NORBT,NORBT,1,LUPRI)
                     ENDIF
                  END IF
               END IF
C
C if both C and D are inactive we also have
C
C        - 3*(PC|C^D) -> FCSO(P,D)
C        + 3*(PC|C^D) -> FCSO(D,P)
C
               IF ( IDPSYM.EQ.KSYMOP .AND. ITYPCD.EQ.JTININ ) THEN
                  IF (IPRHSO.GT.IPRFSO) THEN
                     WRITE(LUPRI,'(A,I3,A,I3,I5)')' PD contribution ',
     *                  IOFFP+1,':',IOFFP+NDIMP,IDI
                     CALL OUTPUT(H2(IOFFP+1,ICI),1,NORBP,1,1,
     *                                         NORBT,NORBT,1,LUPRI)
                  ENDIF
                  CALL DAXPY(NORBP,-FAC,H2(IOFFP+1,ICI),1,
     *                                 FCSO(IOFFP+1,IDI),1)
                  CALL DAXPY(NORBP,FAC,H2(IOFFP+1,ICI),1,
     *                                 FCSO(IDI,IOFFP+1),NORBT)
               ENDIF
            ENDIF
C
         IF (LORB) THEN
C
C Active Fock matrix
C
C Direct terms:
C
C  FVSO(P,Q) = FVSO(P,Q)+ SUM(X,Y) (XY|P^Q)*D(XY)
C
C
C here:
C        + SUM(X,Y) (XY|C^D)*D(XY) -> FVSO(C,D)
C        - SUM(X,Y) (XY|C^D)*D(XY) -> FVSO(D,C)
C
C where the sum is taken over diagonal symmetry blocks (X,Y)
C
            FAC = DISFAC*COUFAC
            IXSYM = ISYM
            IASHX = IASH(IXSYM)
            NASHX = NASH(IXSYM)
            IOFFXA = IORB(IXSYM) + NISH(IXSYM)
            IF (KSYMOP.EQ.ICDSYM) THEN
               IF (IPRHSO.GT.IPRFSO) THEN
                  WRITE(LUPRI,'(/A)')' Active direct terms'
               END IF
               DO 20 IX=1,NASHX
                  WRK(IX) = DDOT(NASHX,H2(IOFFXA+1,IOFFXA+IX),1,
     *                           UDV(IASHX+1,IASHX+IX),1)
 20            CONTINUE
               XYSUM = DSUM(NASHX,WRK(1),1)
               FVSO(ICI,IDI) = FVSO(ICI,IDI) + FAC*XYSUM
               FVSO(IDI,ICI) = FVSO(IDI,ICI) - FAC*XYSUM
            END IF
C
C Exchange terms: rearranged with C position active
C
C  FVSO(P,Q) = FVSO(P,Q) - SUM(X,Y) 3/2*(PX|Y^Q)*D(X,Y)
C                        + SUM(X,Y) 3/2*(QX|Y^P)*D(X,Y)
C
C
C here:
C       - SUM(X) 3/2*(PX|C^D)*D(X,C) -> FVSO(P,D)
C       + SUM(X) 3/2*(PX|C^D)*D(X,C) -> FVSO(D,P)
C
C for active-active and active-inactive distributions
C
            FAC=D1P5*DISFAC
            IXSYM = MULD2H(IPSYM,KCDSYM)
            IF (IPRHSO.GT.IPRFSO) THEN
               IF (IXSYM.EQ.ICSYM .AND.
     *                   (ITYPCD.EQ.JTACIN .OR. ITYPCD.EQ.JTACAC)
     *                          .OR.
     *             IXSYM.EQ.IDSYM .AND.
     *                   (ITYPCD.EQ.JTACAC .OR.
     *                          NOCCP.GT.0 .AND. ITYPCD.EQ.JTSEAC)) THEN
                  WRITE(LUPRI,'(/A)')' Active exchange terms'
               END IF
            END IF
            IF (IXSYM.EQ.ICSYM) THEN
               IOFFXA = IORB(IXSYM) + NISH(IXSYM)
               NASHX = NASH(IXSYM)
               IASHX = IASH(IXSYM)
               IF (ITYPCD.EQ.JTACIN .OR. ITYPCD.EQ.JTACAC) THEN
                  CALL DGEMM('N','N',NORBP,1,NASHX,1.0D0,
     &                       H2(IOFFP+1,IOFFXA+1),NORBT,
     &                       UDV(IASHX+1,NCIW),NASHT,0.0D0,WRK(1),NORBP)
                  CALL DAXPY(NORBP,-FAC,WRK(1),1,
     *                              FVSO(IOFFP+1,IDI),1)
                  CALL DAXPY(NORBP,FAC,WRK(1),1,
     *                             FVSO(IDI,IOFFP+1),NORBT)
                  IF (IPRHSO.GT.IPRFSO) THEN
                     WRITE(LUPRI,'(A,I3,A,I3,I5)')' PD contribution ',
     *                  IOFFP+1,':',IOFFP+NORBP,IDI
                     CALL OUTPUT(WRK(1),1,NORBP,1,1,NORBP,1,1,LUPRI)
                  ENDIF
               END IF
            END IF
C
C Exchange terms: rearranged with D position active
C
C  FVSO(P,Q) = FVSO(P,Q) + SUM(X,Y) 3/2*(PX|Q^Y)*D(X,Y)
C                        - SUM(X,Y) 3/2*(QX|P^Y)*D(X,Y)
C
C
C here:
C        SUM(X) 3/2*(PX|C^D)*D(X,D) -> FVSO(P,C)
C      - SUM(X) 3/2*(PX|C^D)*D(X,D) -> FVSO(C,P)
C
C such that (CP) is at most secondary-active
C
C for active-active and secondary-active distributions
C
            IF (IXSYM.EQ.IDSYM) THEN
               IOFFXA = IORB(IXSYM) + NISH(IXSYM)
               NASHX = NASH(IXSYM)
               IASHX = IASH(IXSYM)
               IF (ITYPCD.EQ.JTACAC .OR. ITYPCD.EQ.JTSEAC) THEN
                  IF (ITYPCD.EQ.JTSEAC) THEN
                     NDIMP = NOCCP
                  ELSE
                     NDIMP = NORBP
                  END IF
                  IF (NDIMP .GT. 0) THEN
                     CALL DGEMM('N','N',NDIMP,1,NASHX,1.0D0,
     &                          H2(IOFFP+1,IOFFXA+1),NORBT,
     &                          UDV(IASHX+1,NDIW),NASHT,0.0D0,
     &                          WRK(1),NDIMP)
                     CALL DAXPY(NDIMP,FAC,WRK(1),1,
     *                          FVSO(IOFFP+1,ICI),1)
                     CALL DAXPY(NDIMP,-FAC,WRK(1),1,
     *                          FVSO(ICI,IOFFP+1),NORBT)
                     IF (IPRHSO.GT.IPRFSO) THEN
                       WRITE(LUPRI,'(A,I3,A,I3,I5)')' PC contribution ',
     *                       IOFFP+1,':',IOFFP+NDIMP,ICI
                       CALL OUTPUT(WRK(1),NDIMP,1,1,1,NDIMP,1,1,LUPRI)
                     ENDIF
                  END IF
               END IF
            END IF
C
         END IF
C        (RSPCI)
200      CONTINUE
9999  CALL QEXIT('FSOMU')
      RETURN
      END
C  /* Deck qsomu */
      SUBROUTINE QSOMU(ICI,IDI,QASO,QBSO,
     *                  H2,PVX,PV12,PV21,
     *                  WRK,LWRK)
C
C Olav Vahtras
C Apr 11, 1990
C
C Purpose:
C  Calculate all contributions to QA and QB spin-orbit matrices
C  from Mulliken (**|C^D) integral distributions
C
C  In general:
C
C  QBSO(P,Q) = SUM(X,Y,W) (PW|X^Y)*( 2*PV(++)(Q,W,X,Y) + PV(--)(Q,W,X,Y) )
C            + SUM(X,Y,W) (XY|P^W)*( 2*PV(--)(X,Y,Q,W) + PV(++)(X,Y,Q,W) )
C
C  QASO(P,Q) = SUM(X,Y,W) (WP|X^Y)*( 2*PV(++)(W,Q,X,Y) + PV(--)(W,Q,X,Y) )
C            + SUM(X,Y,W) (XY|W^P)*( 2*PV(--)(X,Y,W,Q) + PV(++)(X,Y,W,Q) )
C
C   where PV(++)(P,Q,R,S) = <0| e(+,+)(pqrs) |0>
C    and  PV(--)(P,Q,R,S) = <0| e(-,-)(pqrs) |0>
C
#include "implicit.h"
C
#include "maxash.h"
#include "maxorb.h"
#include "priunit.h"
#include "inforb.h"
#include "infind.h"
#include "wrkrsp.h"
#include "infrsp.h"
#include "infpri.h"
#include "orbtypdef.h"
#include "trhso.h"
#include "infhso.h"
C
C
      DIMENSION QASO(NORBT,NASHT),QBSO(NORBT,NASHT)
      DIMENSION PV12(NASHT,NASHT),PV21(NASHT,NASHT)
      DIMENSION H2(NORBT,NORBT),PVX(*),WRK(*)
C
      PARAMETER(D2=2.0D0)
C
      CALL QENTER('QSOMU')
C
C  Local print level
C
      IPRQSO = 20
C
      IF (LWRK.LT.N2ASHX) CALL ERRWRK('QSOMU',N2ASHX,LWRK)
C
C
C Symmetry to type order
C
      ICIW = ISW(ICI)
      IDIW = ISW(IDI)
C
C Orbital type, at least one has to be active to contribute
C
      ITYPC = IOBTYP(ICI)
      ITYPD = IOBTYP(IDI)
      IF (ITYPC.NE.JTACT .AND. ITYPD.NE.JTACT) GO TO 9999
C
C Distribution type and symmetry
C
      ICSYM = ISMO(ICI)
      IDSYM = ISMO(IDI)
      ICDSYM = MULD2H(ICSYM,IDSYM)
      KCDSYM = MULD2H(KSYMOP,ICDSYM)
      ITYPCD=IDBTYP(ITYPC,ITYPD)
C
C     Order within actives for the case that C or D are active
C
      IF (ITYPC .EQ. JTACT) NCIW = ICIW - NISHT
      IF (ITYPD .EQ. JTACT) NDIW = IDIW - NISHT
C
      IF (IPRHSO.GT.IPRQSO) THEN
         WRITE(LUPRI,'(/A)') ' ------ Output from QSOMU ------'
         WRITE(LUPRI,'(/A,2I5,5X,2A)')' Distribution CD',ICI,IDI,
     *                          COBTYP(ITYPC),COBTYP(ITYPD)
         WRITE(LUPRI,'(A,2I5)')       ' Symmetry       ',ICSYM,IDSYM
      ENDIF
C
C
            IPP = 1
            IMM = N2ASHX*N2ASHX+1
C
C Both C and D are active:
C
      IF (ITYPCD.EQ.JTACAC) THEN
C
C Get (C,D) density distributions in the form  2PV(++)+PV(--)
C
            NCDOFF = (NCIW-1 + (NDIW-1)*NASHT)*N2ASHX
            CALL DCOPY(N2ASHX,PVX(NCDOFF+IMM),1,PV12,1)
            CALL DAXPY(N2ASHX,D2,PVX(NCDOFF+IPP),1,PV12,1)
            IF (IPRHSO.GT.IPRQSO) THEN
               WRITE(LUPRI,'(/A)')' Active diagonal terms'
               WRITE(LUPRI,'(/A,2I5)')
     *        ' Total CD density distribution' ,NCIW,NDIW
               CALL OUTPUT(PV12,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI)
            END IF
C
C Get (D,C) density distributions
C
            NDCOFF = (NDIW-1 + (NCIW-1)*NASHT)*N2ASHX
            CALL DCOPY(N2ASHX,PVX(NDCOFF+IMM),1,PV21,1)
            CALL DAXPY(N2ASHX,D2,PVX(NDCOFF+IPP),1,PV21,1)
            IF (IPRHSO.GT.IPRQSO) THEN
               WRITE(LUPRI,'(/A,2I5)')
     *       ' Total DC density distribution' ,NDIW,NCIW
               CALL OUTPUT(PV21,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI)
            END IF
C
C Add contibutions to QASO and QBSO from (C,D) and (D,C) density
C distributions
C
C
               IF (IPRHSO.GT.IPRQSO) THEN
                  WRITE(LUPRI,'(/A)') ' Loop over symmetry blocks '
               END IF
               DO 100 IPSYM = 1,NSYM
                  IWSYM = MULD2H(IPSYM,KCDSYM)
                  IQSYM = MULD2H(IWSYM,ICDSYM)
                  NORBP = NORB(IPSYM)
                  NASHP = NASH(IPSYM)
                  NASHQ = NASH(IQSYM)
                  NASHW = NASH(IWSYM)
                  IF (IPRHSO.GT.IPRQSO) THEN
                     WRITE(LUPRI,'(/A,I5)')' IPSYM',IPSYM
                  END IF
                  IF ( NORBP.GT.0 .AND. NASHW.GT.0 .AND. NASHQ.GT.0 )
     *                                                           THEN
                     IOFFP = IORB(IPSYM)
                     IOFFW = IORB(IWSYM)
                     IOFFPA = IOFFP + NISH(IPSYM)
                     IOFFWA = IOFFW + NISH(IWSYM)
                     IASHP = IASH(IPSYM)
                     IASHQ = IASH(IQSYM)
                     IASHW = IASH(IWSYM)
                     IF (IPRHSO.GT.IPRQSO) THEN
                        WRITE(LUPRI,'(/A,2I5)')
     *                  ' Integral symmetry block PW',IPSYM,IWSYM
                        CALL OUTPUT(H2,IOFFP+1,IOFFP+NORBP,
     *                                 IOFFWA+1,IOFFWA+NASHW,
     *                                 NORBT,NORBT,1,LUPRI)
                        WRITE(LUPRI,'(/A,2I5)')
     *                  ' Density symmetry block QW',IQSYM,IWSYM
                        CALL OUTPUT(PV12,IASHQ+1,IASHQ+NASHQ,
     *                                   IASHW+1,IASHW+NASHW,
     *                                   NASHT,NASHT,1,LUPRI)
                     END IF
C
C QBSO(P,Q) += SUM(W) (PW|C^D)*(2*PV(++)(Q,W,C,D) + PV(--)(Q,W,C,D))
C
C
               CALL DGEMM('N','T',NORBP,NASHQ,NASHW,1.D0,
     &                    H2(IOFFP+1,IOFFWA+1),NORBT,
     &                    PV12(IASHQ+1,IASHW+1),NASHT,1.D0,
     &                    QBSO(IOFFP+1,IASHQ+1),NORBT)
C
C QBSO(P,Q) += SUM(W) (PW|D^C)*(2*PV(++)(Q,W,D,C) + PV(--)(Q,W,D,C))
C        (  -     ...     C^D   )
C
               CALL DGEMM('N','T',NORBP,NASHQ,NASHW,-1.D0,
     &                    H2(IOFFP+1,IOFFWA+1),NORBT,
     &                    PV21(IASHQ+1,IASHW+1),NASHT,1.D0,
     &                    QBSO(IOFFP+1,IASHQ+1),NORBT)
C
C QASO(P,Q) += SUM(W) (WP|C^D)*(2*PV(++)(W,Q,C,D) + PV(--)(W,Q,C,D))
C
               CALL DGEMM('T','N',NORBP,NASHQ,NASHW,1.D0,
     &                    H2(IOFFWA+1,IOFFP+1),NORBT,
     &                    PV12(IASHW+1,IASHQ+1),NASHT,1.D0,
     &                    QASO(IOFFP+1,IASHQ+1),NORBT)
C
C QASO(P,Q) += SUM(W) (WP|D^C)*(2*PV(++)(W,Q,D,C) + PV(--)(W,Q,D,C))
C        (  -     ...     C^D   )
C
               CALL DGEMM('T','N',NORBP,NASHQ,NASHW,-1.D0,
     &                    H2(IOFFWA+1,IOFFP+1),NORBT,
     &                    PV21(IASHW+1,IASHQ+1),NASHT,1.D0,
     &                    QASO(IOFFP+1,IASHQ+1),NORBT)
C
                  END IF
 100           CONTINUE
      END IF
C
C Following integral distributions contribute if either
C C or D is active.
C
C
C QBSO(C,Q) += SUM(X,Y) (XY|C^D)*(2*PV(--)(X,Y,Q,D)+PV(++)(X,Y,Q,D))
C
      IF (ITYPD.EQ.JTACT) THEN
         IQSYM = MULD2H(KSYMOP,ICSYM)
         NASHQ = NASH(IQSYM)
         IASHQ = IASH(IQSYM)
         IF (NASHQ.GT.0) THEN
            IF (IPRHSO.GT.IPRQSO) THEN
               WRITE(LUPRI,'(/A,I5)')
     *         ' Loop over actives in symmetry',IQSYM
            END IF
            DO 200 IQ=1,NASHQ
C
C For each Q read a new (Q,D) density distribution in the form
C 2PV(--) + PV(++)
C
               IF (IPRHSO.GT.IPRQSO) THEN
                  WRITE(LUPRI,'(/A,2I5)')' IQ',IASHQ+IQ
               END IF
               NQDOFF = (IASHQ+IQ-1 + (NDIW-1)*NASHT)*N2ASHX
               CALL DCOPY(N2ASHX,PVX(NQDOFF+IPP),1,PV12,1)
               CALL DAXPY(N2ASHX,D2,PVX(NQDOFF+IMM),1,PV12,1)
               IF (IPRHSO.GT.IPRQSO) THEN
                  WRITE(LUPRI,'(/A,2I5)')
     *           ' Total QD density distribution', IASHQ+IQ,NDIW
                  CALL OUTPUT(PV12,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI)
               END IF
C
C
                  DO 20 IYSYM = 1,NSYM
                     NASHY = NASH(IYSYM)
                     IASHY = IASH(IYSYM)
                     IOFFYA = IORB(IYSYM) + NISH(IYSYM)
                     IXSYM = MULD2H(IYSYM,KCDSYM)
                     NASHX = NASH(IXSYM)
                     IASHX = IASH(IXSYM)
                     IOFFXA = IORB(IXSYM) + NISH(IXSYM)
                     IF (NASHX.GT.0 .AND .NASHY .GT.0) THEN
C
C QBSO(C,Q) += SUM(X,Y) (XY|C^D)*(2*PV(--)(X,Y,Q,D)+PV(++)(X,Y,Q,D))
C
C
C QASO(C,Q) += SUM(X,Y) (XY|D^C)*(2*PV(--)(X,Y,D,Q)+PV(++)(X,Y,D,Q))
C        (  -       ...     C^D   )           (Q,D)           (Q,D)
C
                     IF (IPRHSO.GT.IPRQSO) THEN
                        WRITE(LUPRI,'(/A,2I5)')
     *                  ' Density symmetry block XY',IXSYM,IYSYM
                        CALL OUTPUT(PV12,IASHX+1,IASHX+NASHX,
     *                                   IASHY+1,IASHY+NASHY,
     *                                   NASHT,NASHT,1,LUPRI)
                     END IF
                        DO 21 IY=1,NASHY
                          WRK(IY) = DDOT(NASHX,H2(IOFFXA+1,IOFFYA+IY),1,
     *                                         PV12(IASHX+1,IASHY+IY),1)
 21                     CONTINUE
                        XYSUM = DSUM(NASHY,WRK,1)
                        QBSO(ICI,IASHQ+IQ)
     *                 = QBSO(ICI,IASHQ+IQ) + XYSUM
                        QASO(ICI,IASHQ+IQ)
     *                 = QASO(ICI,IASHQ+IQ) - XYSUM
                     END IF
 20               CONTINUE
 200        CONTINUE
         END IF
      ENDIF
C
C QBSO(D,Q) += SUM(X,Y) (XY|D^C)*(2*PV(--)(X,Y,Q,C)+PV(++)(X,Y,Q,C))
C
      IF (ITYPC.EQ.JTACT) THEN
         IQSYM = MULD2H(KSYMOP,IDSYM)
         NASHQ = NASH(IQSYM)
         IASHQ = IASH(IQSYM)
         IF (NASHQ.GT.0) THEN
            IF (IPRHSO.GT.IPRQSO) THEN
               WRITE(LUPRI,'(/A,I5)')
     *         ' Loop over actives in symmetry',IQSYM
            END IF
            DO 300 IQ = 1,NASHQ
               IF (IPRHSO.GT.IPRQSO) THEN
                  WRITE(LUPRI,'(/A,2I5)')' IQ',IASHQ+IQ
               END IF
C
C For each Q read a new (Q,C) density distribution in the form
C 2PV(--) + PV(++)
C
               NQCOFF = (IASHQ+IQ-1 + (NCIW-1)*NASHT)*N2ASHX
               CALL DCOPY(N2ASHX,PVX(NQCOFF+IPP),1,PV21,1)
               CALL DAXPY(N2ASHX,D2,PVX(NQCOFF+IMM),1,PV21,1)
               IF (IPRHSO.GT.IPRQSO) THEN
                  WRITE(LUPRI,'(/A,2I5)')
     *           ' Total QC density distribution', IASHQ+IQ,NCIW
                  CALL OUTPUT(PV21,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI)
               END IF
C
C
                  DO 30 IYSYM = 1,NSYM
                     NASHY = NASH(IYSYM)
                     IASHY = IASH(IYSYM)
                     IOFFYA = IORB(IYSYM) + NISH(IYSYM)
                     IXSYM = MULD2H(IYSYM,KCDSYM)
                     NASHX = NASH(IXSYM)
                     IASHX = IASH(IXSYM)
                     IOFFXA = IORB(IXSYM) + NISH(IXSYM)
                     IF (NASHX.GT.0 .AND. NASHY.GT.0) THEN
C
C QBSO(D,Q) += SUM(X,Y) (XY|D^C)*(2*PV(--)(X,Y,Q,C)+PV(++)(X,Y,Q,C))
C        (  -       ...     C^D   )
C
C
C QASO(D,Q) += SUM(X,Y) (XY|C^D)*(2*PV(--)(X,Y,C,Q)+PV(++)(X,Y,C,Q))
C                                             (Q,C)           (Q,C)
                     IF (IPRHSO.GT.IPRQSO) THEN
                        WRITE(LUPRI,'(/A,2I5)')
     *                  ' Density symmetry block XY',IXSYM,IYSYM
                        CALL OUTPUT(PV21,IASHX+1,IASHX+NASHX,
     *                                   IASHY+1,IASHY+NASHY,
     *                                   NASHT,NASHT,1,LUPRI)
                     END IF
                        DO 31 IY=1,NASHY
                          WRK(IY) = DDOT(NASHX,H2(IOFFXA+1,IOFFYA+IY),1,
     *                                         PV21(IASHX+1,IASHY+IY),1)
 31                     CONTINUE
                        XYSUM = DSUM(NASHY,WRK,1)
                        QBSO(IDI,IASHQ+IQ)
     *                 = QBSO(IDI,IASHQ+IQ) - XYSUM
                        QASO(IDI,IASHQ+IQ)
     *                 = QASO(IDI,IASHQ+IQ) + XYSUM
                     END IF
 30               CONTINUE
C
 300        CONTINUE
         END IF
      END IF
C
 9999 CALL QEXIT('QSOMU')
      RETURN
      END
C  /* Deck hsosig */
      SUBROUTINE HSOSIG(FC,H2AC,
     *                  GP,XINDX,WRK,LWRK)
C
C  Calculate the configuration part of the spin-orbit property vector
C
C                ( <J,HSO,0> )
C      GP(J) =   (           )
C                (-<0,HSO,J> )
C
C
C H2AC contain active two-electron spin-orbit integrals in doubly
C triangularly packed form
C
C Olav Vahtras
C Sep 14, 1990
#include "implicit.h"
C
      DIMENSION FC(NORBT,NORBT)
      DIMENSION H2AC(NNASHX,NNASHX)
      DIMENSION GP(KZYVAR),XINDX(*),WRK(*)
C
#include "maxorb.h"
#include "maxash.h"
#include "priunit.h"
#include "wrkrsp.h"
#include "infrsp.h"
#include "infopt.h"
#include "inforb.h"
#include "infind.h"
#include "infpri.h"
#include "infdim.h"
#include "cbgetdis.h"
#include "trhso.h"
#include "infhso.h"
C
C
      PARAMETER ( DM1 = -1.0D0 )
C
      CALL QENTER('HSOSIG')
C
C  Local print level
C
      IPRLOC = 20
      IF (IPRHSO.GT.IPRLOC) THEN
         WRITE(LUPRI,'(/A)') ' ------ Output from HSOSIG ------'
      END IF
C
C ALLOCATE WORK SPACE
C
      IF (IREFSY .EQ. KSYMST) THEN
         NDREF = KZCONF
C        ... if KZCONF .ne. NCREF, we need CREF in determinants
      ELSE
         NDREF = NCREF
      END IF
      KFREE = 1
      LFREE = LWRK
      CALL MEMGET2('REAL','FCAC' ,KFCAC,N2ASHX,WRK,KFREE,LFREE)
      CALL MEMGET2('REAL','REFCO',KREFCO,NDREF,WRK,KFREE,LFREE)
      CALL MEMGET2('REAL','HSQSQ',KHSQSQ,N2ASHX*N2ASHX,WRK,KFREE,LFREE)
      CALL DZERO(WRK(KFCAC),N2ASHX)
      CALL DZERO(WRK(KREFCO),NDREF)
      CALL DZERO(WRK(KHSQSQ),N2ASHX*N2ASHX)
      IF (TESTZY) THEN
         CALL MEMGET2('REAL','HSQTR',KHSQTR,N2ASHX*N2ASHX,
     &      WRK,KFREE,LFREE)
         CALL DZERO(WRK(KHSQTR),N2ASHX*N2ASHX)
      END IF
      IADINT = -1
C
C
C Obtain square packed combined two-electron integrals
C
      CALL H2ACSO(H2AC,WRK(KHSQSQ))
C
      ISPIN1 = 0
      ISPIN2 = 1
C
      CALL GETREF(WRK(KREFCO),NDREF)
C
C CREATE Z PART  <0,HSO,J> OF SPIN-ORBIT GP VECTOR
C
C        Get FCAC matrix for Z sigma vector
C        (note: CISIGD requires UFCAC(I,J) = FCXAC(J,I))
C
         DO 220 IW = 1,NASHT
            IX = ISX(NISHT+IW)
            DO 230 JW = 1,NASHT
               JX  = ISX(NISHT+JW)
               IJW = (IW-1) * NASHT + JW
               WRK(KFCAC-1+IJW) = FC(IX,JX)
 230        CONTINUE
 220     CONTINUE
C
         IF (IPRHSO.GT.IPRLOC) THEN
            WRITE(LUPRI,'(/A)')' Inactive Fock matrix'
            CALL OUTPUT(FC,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
            WRITE(LUPRI,'(/A)')' Active part of inactive Fock matrix'
            CALL OUTPUT(WRK(KFCAC),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI)
         END IF
         DISTYP = 6
C
         KSYMST = MULD2H(IREFSY,KSYMSO)
         CALL CISIGD(IREFSY, KSYMST, NDREF, KZCONF, WRK(KREFCO),GP,
     *               WRK(KFCAC), WRK(KHSQSQ), .FALSE., .FALSE.,
     *               XINDX, ISPIN1, ISPIN2, WRK(KFREE), LFREE)
C        CALL RSPSIG(ICSYM,IHCSYM,NCDET,NHCDET,C,HC,UFCAC,H2AC,IFLAG,
C    *               NOH2,WORK,KFREE,LFREE)
C
C        IF ((NDREF.NE.NCDET).OR.(KZCONF.NE.NHCDET)) THEN
C           WRITE(LUPRI,'(/2(A,I5),/3(A,I5))')
C    *      ' NUMBER OF REFERENCE DETERMINANTS ,NDREF:',NDREF,
C    *      ' CALCULATED NUMBER ,NCDET:',NCDET,
C    *      ' NUMBER OF DETERMINANTS FOR SYMMETRY',KSYMOP,
C    *      '  IS:',KZCONF,'  CALCULATED NUMBER, NHCDET:',NHCDET
C        CALL QUIT(' H2XSIG, INCORRECT CALCULATION OF DETERMINANTS')
C        END IF
         IF (IPRHSO.GT.IPRLOC) THEN
            WRITE(LUPRI,'(/A)')
     *      ' Configuration part of spin-orbit property vector: Z part'
            CALL OUTPUT(GP,1,KZCONF,1,1,KZCONF,1,1,LUPRI)
         END IF
C
C CREATE Y PART  <0,HSO,J> OF SPIN-ORBIT GP VECTOR
C
C Normal procedure: copy the Z part
C Test procedure: do it explicitely
C
C        Get transposed FCXAC matrix for Y sigma vector
C        (note: CISIGD requires UFCAC(I,J) = FCXAC(J,I))
C
         IF (TESTZY) THEN
C
C           ISPIN1 = 0
C           ISPIN2 = 1
C
            DO 120 IW = 1,NASHT
               IX = ISX(NISHT+IW)
               DO 130 JW = 1,NASHT
                  JX = ISX(NISHT+JW)
                  IJW = (IW-1) * NASHT + JW
                  WRK(KFCAC-1+IJW) = FC(JX,IX)
 130           CONTINUE
 120        CONTINUE
            IF (IPRHSO.GT.IPRLOC) THEN
               WRITE(LUPRI,'(/A)')' Inactive Fock matrix'
               CALL OUTPUT(FC,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
               WRITE(LUPRI,'(/A)')' Active part of inactive Fock matrix'
               CALL OUTPUT(WRK(KFCAC),1,NASHT,1,NASHT,NASHT,NASHT,1,
     &                    LUPRI)
            END IF
C
C Get transposed integrals
C
C           CALL MTRSP(N2ASHX,N2ASHX,WRK(KHSQSQ),N2ASHX,
C    *                               WRK(KHSQTR),N2ASHX)
            CALL DAXPY(N2ASHX*N2ASHX,DM1,WRK(KHSQSQ),1,WRK(KHSQTR),1)
            DISTYP = 6
            CALL CISIGD(IREFSY, KSYMST, NDREF, KZCONF, WRK(KREFCO),
     *               GP(KZVAR+1),
     *               WRK(KFCAC), WRK(KHSQTR), .FALSE., .FALSE.,
     *               XINDX, ISPIN1, ISPIN2, WRK(KFREE), LFREE)
            CALL DSCAL(KZCONF,DM1,GP(1+KZVAR),1)
            IF (IPRHSO.GT.IPRLOC) THEN
               WRITE(LUPRI,'(/A)')
     &      ' Configuration part of spin-orbit property vector: Y part'
               CALL OUTPUT(GP(1+KZVAR),1,KZCONF,1,1,KZCONF,1,1,LUPRI)
            END IF
C
         ELSE
            CALL DCOPY(KZCONF,GP,1,GP(1+KZVAR),1)
         END IF
 100  CONTINUE
      CALL MEMREL('HSOSIG',WRK,1,1,KFREE,LFREE)
      CALL QEXIT('HSOSIG')
      RETURN
      END
C  /* Deck h2acso */
      SUBROUTINE H2ACSO(H2ACPK,H2SOAC)
C
C 20-Jun-1990 hjaaj
C
#include "implicit.h"
      DIMENSION H2ACPK(NNASHX,NNASHX), H2SOAC(N2ASHX,NASHT,NASHT)
C
      PARAMETER ( DM1 = -1.0D0 )
C
C Used from common blocks:
C  INFORB : NASHT,NNASHX,N2ASHX
C  INFIND : IROW()
C
#include "maxash.h"
#include "maxorb.h"
#include "inforb.h"
#include "infind.h"
C
      CALL QENTER('H2ACSO')
C
C     1. Unpack H2ACPK into H2SOAC
C
      DO 140 K = 1,NASHT
         DO 130 L = 1,K-1
            KL = IROW(K) + L
            CALL DSPTSI(NASHT,H2ACPK(1,KL),H2SOAC(1,K,L))
            CALL DAXPY(N2ASHX,DM1,H2SOAC(1,K,L),1,H2SOAC(1,L,K),1)
  130    CONTINUE
         CALL DZERO(H2SOAC(1,K,K),N2ASHX)
  140 CONTINUE
      CALL QEXIT('H2ACSO')
      RETURN
      END
      SUBROUTINE TRANSFORM_HSO(WORD,CMO,WRK,KWRK1,LWRK1)
#include "implicit.h"
      CHARACTER*8 WORD
      DIMENSION  WRK(*), CMO(*)
#include "dummy.h"
#include "rspprp.h"
#include "infhso.h"
#include "infhyp.h"
#include "infpp.h"
#include "infsmo.h"
#include "inftap.h"
#include "priunit.h"
#include "trhso.h"
#include "wrkrsp.h"
      LOGICAL FOPEN
      LOGICAL MOEXIS, SO2TRA
C
C Set spin-orbit component ILXYZ for TRHSO
C
      CALL QENTER('TRANSFORM_HSO')
      IF (WORD(1:1).EQ.'X') THEN
         ILXYZ = 1
      ELSE IF (WORD(1:1).EQ.'Y') THEN
         ILXYZ = 2
      ELSE IF (WORD(1:1).EQ.'Z') THEN
         ILXYZ = 3
      ELSE
         CALL QUIT('Wrong property in TRANSFORM_HSO, WORD = '//WORD)
      END IF
C
C Transform spin-orbit two-electron integrals, we need sec-occ in linear
C response and sec-sec in quadratic response
C
      IF (HYPCAL.OR.SOMOM.OR.EXMOM) THEN
         ITRLSO = 2
      ELSE
         ITRLSO = 1
      END IF
      KSYMSO = KSYMOP
C defined in RSPSET (920922-ov)
      SO2TRA = .TRUE.
C
C If .OLDTRA has been specified check that MO2SOINT exists
C If it does not exist transform AO2SOINT as usual
C
      IF (OLDTRA) THEN
         INQUIRE(FILE='MOHSOINT',EXIST=MOEXIS)
         IF (MOEXIS) THEN
            SO2TRA = .FALSE.
         ELSE
            WRITE(LUPRI,'(/3A/A)') ' WARNING: Expected transformed ',
     *                       ' two-electron spin-orbit integrals not',
     *                       ' found.' ,
     *                       ' - file MOHSOINT does not exist'
            NWARN = NWARN + 1
         END IF
      END IF
      IF (SO2TRA) THEN
         IF (IPRHSO.GT.0) THEN
            CALL TIMER('START ',TRASTA,TRATIM)
         END IF
         CALL GPOPEN(LUAHSO,'AO2SOINT','OLD',' ','UNFORMATTED',IDUMMY,
     &               .FALSE.)
         CALL TRAHSO(ITRLSO,CMO,WRK(KWRK1),LWRK1)
         CALL GPCLOSE(LUAHSO,'KEEP')
         IF (IPRHSO.GT.0) CALL TIMER('TRAHSO',TRASTA,TRATIM)
      ELSE
         INQUIRE(FILE='MOHSOINT',OPENED=FOPEN)
         IF (.NOT.FOPEN) THEN
            CALL DAOPEN(LUMHSO,'MOHSOINT')
         END IF
      END IF
      CALL QEXIT('TRANSFORM_HSO')
      END
      SUBROUTINE HSOFCK(WORD,FC,FV,QA,QB,H2AC,UDV,PV,CMO,LORB,LCON,
     &   WRK,KWRK1,LWRK1)
#include "implicit.h"
      CHARACTER*8 WORD
      DIMENSION FC(*), FV(*), QA(*), QB(*), H2AC(*), UDV(*), PV(*),
     &   CMO(*), WRK(*)
      LOGICAL LORB, LCON
#include "infrsp.h"
#include "inforb.h"
      CALL QENTER('HSOFCK')
      IF (TDHF) THEN
         IF (NASHT.GT.0) THEN
            CALL DZERO(QA,NORBT*NASHT)
            CALL DZERO(QB,NORBT*NASHT)
         END IF
         CALL HSOFCKAO(WORD,CMO,UDV,FC,FV,WRK,KWRK1,LWRK1)
      ELSE
         CALL HSOFCKMO(FC,FV,QA,QB,H2AC,UDV,PV,LORB,LCON,
     &      WRK,KWRK1,LWRK1)
      END IF
      CALL QEXIT('HSOFCK')
      END 
C
C /* Deck hsofckao */
C
      SUBROUTINE HSOFCKAO(WORD,CMO,UDV,FC,FV,WRK,KFREE,LFREE)
      IMPLICIT NONE
      CHARACTER*(*) WORD
      REAL*8    CMO(*), UDV(*), FC(*), FV(*), WRK(*)
      INTEGER KFREE,LFREE
C
C Build HSO inactive fockmatrix from AO integrals
C    If DIRFCK=.TRUE. call TWOINT
C  else read from AO2SOINT
C
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "dummy.h"
#include "infinp.h"
#include "wrkrsp.h"
#include "inforb.h"
#include "aovec.h"
#include "infrsp.h"
#include "dftcom.h"
#include "symmet.h"
      REAL*8 D1
      PARAMETER (D1=1.0D0)
      INTEGER NDMAT,KDMAO,KFMAO,KDMSO,I
      INTEGER KDC,KDV,KFC,KFV,KFMO,KFAO
      INTEGER KJSTRS , KNPRIM , KNCONT , KIORBS , KJORBS , NPAO , KKORBS
      INTEGER I2TYP
      INTEGER IPRTWO,IPRNTA,IPRNTB,IPRNTC,IPRNTD
      INTEGER IRNTYP , NUMDIS, IFCTYP(2), ISYMDM(2)
      REAL*8  TIMEND,TIMSTR,SO1WAL,SO1CPU
      REAL*8  SO2WAL,SO2CPU,SOCPU,SOWAL
      REAL*8  HFXSAV
      LOGICAL TKTIME,RTNTWO
      LOGICAL TRIPLET(2)
      CALL QENTER('HSOFCKAO')

      NDMAT=1
      IF (NASHT.GT.0) NDMAT = 2
      CALL MEMGET2('REAL','FC',KFC,NDMAT*N2BASX,WRK,KFREE,LFREE)
      CALL MEMGET2('REAL','DC',KDC,NDMAT*N2BASX,WRK,KFREE,LFREE)
      KFV = KFC + N2BASX
      KDV = KDC + N2BASX
C
C Full density in 1, spin (=active) density in 2, swap for singlet
C
      CALL GTDMSO(UDV,CMO,WRK(KDC),WRK(KDV),WRK(KFREE))
      IF (NASHT.GT.0) THEN
         CALL DAXPY(N2BASX,D1,WRK(KDV),1,WRK(KDC),1)
         IF (.NOT.TRPLET) THEN
            CALL DSWAP(N2BASX,WRK(KDC),1,WRK(KDV),1)
         END IF
      END IF
      IF (HSODEBUG) THEN
         CALL MAOPRI(WRK(KDC),'HSOFCKAO:Input  D1')
         IF (NASHT.GT.0) THEN
            CALL MAOPRI(WRK(KDC+N2BASX),'HSOFCKAO:Input  D2')
         END IF
      END IF
      IF (DIRFCK) THEN
         IF (WORD(1:1).EQ.'X') IFCTYP(1)=1
         IF (WORD(1:1).EQ.'Y') IFCTYP(1)=2
         IF (WORD(1:1).EQ.'Z') IFCTYP(1)=3
         IF (TRPLET) IFCTYP(1) = -IFCTYP(1)
         IFCTYP(2) = - IFCTYP(1)
C
C        Transform density to AO basis 
C
         CALL DSOTAO(WRK(KDC),WRK(KFC),NBAST,0,IPRRSP)
         IF (NASHT.GT.0) THEN
            CALL DSOTAO(WRK(KDV),WRK(KFV),NBAST,0,IPRRSP)
         END IF
         CALL DCOPY(NDMAT*N2BASX,WRK(KFC),1,WRK(KDC),1)
C
C Setup for TWOINT
C
         NPAO = MXSHEL*MXAOVC
         CALL MEMGET('INTE',KJSTRS,NPAO*2,WRK,KFREE,LFREE)     
         CALL MEMGET('INTE',KNPRIM,NPAO*2,WRK,KFREE,LFREE)     
         CALL MEMGET('INTE',KNCONT,NPAO*2,WRK,KFREE,LFREE)     
         CALL MEMGET('INTE',KIORBS,NPAO  ,WRK,KFREE,LFREE)     
         CALL MEMGET('INTE',KJORBS,NPAO  ,WRK,KFREE,LFREE)     
         CALL MEMGET('INTE',KKORBS,NPAO  ,WRK,KFREE,LFREE)     
         CALL PAOVEC(WRK(KJSTRS),WRK(KNPRIM),WRK(KNCONT),
     &               WRK(KIORBS),WRK(KJORBS),WRK(KKORBS),0,
     &               .FALSE.,IPRRSP)
         CALL MEMREL('HERINT.PAOVEC',WRK,1,KJORBS,KFREE,LFREE)
         CALL TIMER('START ',TIMSTR,TIMEND)
         CALL GETTIM(SO1CPU,SO1WAL)
         I2TYP  = 0
         IRNTYP = -20
         IF (HSODEBUG) THEN
            IPRTWO = 3
         ELSE
            IPRTWO = 0
         END IF
         TKTIME = .FALSE.
         CALL DZERO(WRK(KFC),NDMAT*N2BASX)
C        Always full exchange for spin-orbit integrals:
         HFXSAV=HFXFAC
         HFXFAC=D1
         CALL TWOINT(WRK(KFREE),LFREE,WRK(KFREE),
     &               WRK(KFC),WRK(KDC),NDMAT,
     &               IDUMMY, IFCTYP,
     &               DUMMY,IDUMMY,NUMDIS,1,IRNTYP,0,0,.TRUE.,.TRUE.,
     &               .FALSE.,TKTIME,IPRTWO,IPRNTA,IPRNTB,IPRNTC,
     &               IPRNTD,RTNTWO,IDUMMY,I2TYP,WRK(KJSTRS),
     &               WRK(KNPRIM),WRK(KNCONT),WRK(KIORBS),
     &               IDUMMY,IDUMMY,DUMMY,DUMMY,DUMMY,DUMMY,
     &               .FALSE.,.false.)
         HFXFAC=HFXSAV
         CALL MEMREL('HSOFCKAO.TWOINT',WRK,KFC,KJSTRS,KFREE,LFREE)
         IF (HSODEBUG) THEN
            CALL MAOPRI(WRK(KFC),':TWOINT:Calculated F1')
            IF (NASHT.GT.0) THEN
               CALL MAOPRI(WRK(KFC+N2BASX),':TWOINT:Calculated F2')
            END IF
         END IF
C
C  Transform to SO
C
         ISYMDM(1)=0
         ISYMDM(2)=0
         CALL SKLFCK(WRK(KFC),DUMMY,WRK(KFREE),LFREE,
     &               IPRTWO,.FALSE.,
     &               .FALSE.,.FALSE.,.FALSE.,.TRUE.,IDUMMY,.TRUE.,NDMAT,
     &               ISYMDM,IFCTYP,IDUMMY,.TRUE.)
C
         IF (HSODEBUG) THEN
            CALL MAOPRI(WRK(KFC),':SKLFCK:Calculated F1')
            IF (NASHT.GT.0) THEN
               CALL MAOPRI(WRK(KFC+N2BASX),':SKLFCK:Calculated F2')
            END IF
         END IF
         CALL GETTIM(SO2CPU,SO2WAL)
         SOCPU=SO2CPU-SO1CPU
         SOWAL=SO2WAL-SO1WAL
         CALL TIMER('TWOINT',TIMSTR,TIMEND)
         CALL FLSHFO(LUPRI)
         WRITE(LUPRI,'(2(/A),2(/A,F5.0,A))')
     &      '   Two-electron spin-orbit integrals',
     &   '   =================================',
     &   '   Spin-orbit 2-electron CPU  time ',SOCPU,' seconds',
     &   '   Spin-orbit 2-electron wall time ',SOWAL,' seconds'
         CALL MEMCHK('HSOFCKAO.SKLFCK',WRK,KFC)
      ELSE
C
C Read AO integrals from disk (AO2SOINT)
C
         TRIPLET(1) = TRPLET
         TRIPLET(2) = .NOT.TRPLET
         CALL DZERO(WRK(KFC),NDMAT*N2BASX)
         CALL GET_FSO_AO(WORD,TRIPLET,WRK(KFC),WRK(KDC),NDMAT)
      END IF
      IF (HSODEBUG) THEN
         CALL MAOPRI(WRK(KFC),'HSOFCKAO:Calculated F1')
         IF (NASHT.GT.0) THEN
            CALL MAOPRI(WRK(KFC+N2BASX),'HSOFCKAO:Calculated F2')
         END IF
      END IF
C
C Reorder matrices to match gradient routine
C (fa+fb,fa-fb)/2 = (fc,fo-fc) -> (fo,fc-fo)
C
      IF (NASHT.GT.0) THEN
         CALL DAXPY(N2BASX,D1,WRK(KFV),1,WRK(KFC),1) !now (fo,fo-fc)
C        CALL DSWAP(N2BASX,WRK(KFC),1,WRK(KFV),1)
         CALL DSCAL(N2BASX,-D1,WRK(KFV),1) !now (fo,fc-fo)
      END IF
C
      CALL DZERO(FC,N2ORBX)
      CALL AOTOMO(WRK(KFC),FC,CMO,KSYMOP,WRK(KFREE),LFREE)
      IF (NASHT.GT.0) THEN
         CALL DZERO(FV,N2ORBX)
         CALL AOTOMO(WRK(KFV),FV,CMO,KSYMOP,WRK(KFREE),LFREE)
      END IF
      IF (HSODEBUG) THEN
         CALL MAOPRI(WRK(KFC),'HSOFCKAO:FC(AO)')
         IF (NASHT.GT.0) CALL MAOPRI(WRK(KFV),'HSOFCKAO:FC-FO(AO)')
         CALL MMOPRI(FC,'HSOFCKAO:FO(MO)')
         IF (NASHT.GT.0) CALL MMOPRI(FV,'HSOFCKAO:FC-FO(MO)')
      END IF
      CALL MEMREL('HSOFCKAO',WRK,KFC,KFC,KFREE,LFREE)
      CALL QEXIT('HSOFCKAO')
      END
C  /* Deck hsofckmo  */
      SUBROUTINE HSOFCKMO(FC,FV,QA,QB,H2AC,UDV,PV,LORB,LCON,
     &   WRK,KWRK1,LWRK1)
#include "implicit.h"
      DIMENSION FC(*), FV(*), QA(*), QB(*), H2AC(*), UDV(*), PV(*), 
     &          WRK(*)
      LOGICAL   LORB, LCON
#include "infhso.h"
#include "maxorb.h"
#include "maxash.h"
#include "infind.h"
#include "inforb.h"
#include "priunit.h"
#include "wrkrsp.h"
      INTEGER NEEDSO(-4:6)
      CALL QENTER('HSOFCKMO')
      CALL DZERO(FC,N2ORBX)
      IF (LORB) THEN
         CALL DZERO(FV,N2ORBX)
         CALL DZERO(QA,NORBT*NASHT)
         CALL DZERO(QB,NORBT*NASHT)
      END IF
      IF (LCON) THEN
         CALL DZERO(H2AC,NNASHX*NNASHX)
      END IF
      IF (.NOT.DOSO2) THEN
         CALL QEXIT('HSOFCK')
         RETURN
      END IF
      CALL MEMGET2('REAL','H2'  ,KH2  ,N2ORBX,WRK,KWRK1,LWRK1)
      CALL MEMGET2('REAL','PV12',KPV12,N2ASHX,WRK,KWRK1,LWRK1)
      CALL MEMGET2('REAL','PV21',KPV21,N2ASHX,WRK,KWRK1,LWRK1)
C
C Read distributions
C
      NEEDSO(-4:6) = 0
      NEEDSO( 1:5) = 1
      IDIST = 0
      KFREE = 1
      LFREE = LWRK1
 90   CALL NXTHSO(IC,ID,WRK(KH2),NEEDSO,WRK(KWRK1),KFREE,LFREE,IDIST)
      IF (IDIST.LT.0) GOTO 95
      IF (IC.EQ.ID) GOTO 90
      IF (IPRHSO.GT.20) THEN
         WRITE(LUPRI,'(//A,2I5)')' Integral distribution ',IC,ID
         CALL OUTPUT(WRK(KH2),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
      END IF
C
C   Construct inactive and active spin-orbit Fock matrices
C
      CALL FSOMU(IC,ID,WRK(KH2),FC,FV,UDV,LORB,
     * WRK(KWRK1-1+KFREE),LFREE)
C     CALL FSOMU(ICI1,IDI1,H2,FCSO,FVSO,UDV,LORB,WRK,LWRK)
C
C
C   Construct spin-orbit Q matrices
C
      IF (NASHT.GT.0 .AND. LORB)
     *   CALL QSOMU(IC,ID,QA,QB,WRK(KH2),
     *           PV,WRK(KPV12),WRK(KPV21),
     *           WRK(KWRK1-1+KFREE),LFREE)
C
C Add active-active distributions to H2AC(uv,xy)
C     CALL ADH2AC(H2ACXY,H2XY,IUVSYM)
C
      IF (IOBTYP(IC).EQ.JTACT .AND. IOBTYP(ID).EQ.JTACT .AND.
     *      LCON) THEN
         ICSYM = ISMO(IC)
         IDSYM = ISMO(ID)
         ICDSYM = MULD2H(ICSYM,IDSYM)
         KCDSYM = MULD2H(KSYMOP,ICDSYM)
         NCW = ICH(IC)
         NDW = ICH(ID)
         IF (NCW.GT.NDW) THEN
            NCDW = IROW(NCW) + NDW
            KH2XY = 1 + (NCDW-1)*NNASHX
            CALL ADH2AC(H2AC(KH2XY),WRK(KH2),KCDSYM)
         ELSE
            NCDW = IROW(NDW) + NCW
            KH2XY = 1 + (NCDW-1)*NNASHX
            CALL DAXPY(N2ORBX,-1.0D0,WRK(KH2),1,WRK(KWRK1-1+KFREE),1)
            CALL ADH2AC(H2AC(KH2XY),WRK(KWRK1-1+KFREE),KCDSYM)
         END IF
      END IF
      GOTO 90
 95   CONTINUE
C
      IF (IPRHSO.GT.10) THEN
         WRITE(LUPRI,'(//A)') ' Final Inactive Fock matrix'
         CALL OUTPUT(FC,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
C
         IF (LORB .AND. NASHT.GT.0) THEN
            WRITE(LUPRI,'(//A)') ' Final Active Fock matrix'
            CALL OUTPUT(FV,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
C
            WRITE(LUPRI,'(//A)') ' Spin-orbit QA matrix'
            CALL OUTPUT(QA,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI)
            WRITE(LUPRI,'(//A)') ' Spin-orbit QB matrix'
            CALL OUTPUT(QB,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI)
         END IF
C
         IF (LCON) THEN
            DO 71 K=1,NASHT
               DO 72 L=1,K
                  KL = IROW(K)+L
                  KLOFF = (KL-1)*NNASHX
                  WRITE(LUPRI,'(/A,2I3,A)')
     *         ' Active integral distribution (**', K,L,')'
                  CALL OUTPAK(H2AC(1+KLOFF),NASHT,1,LUPRI)
 72            CONTINUE
 71         CONTINUE
         END IF
      END IF
      CALL MEMREL('HSOFCKMO',WRK,KH2,KH2,KWRK1,LWRK1)
      CALL QEXIT('HSOFCKMO')
      END
C  /* Deck hsoorb */
      SUBROUTINE HSOORB(ONEIND,NSIM,FC,FCX,FVX,QAX,QBX,UDV,
     *                  EVECS,TRPDEN)
C
C WRITTEN 14-FEB 1986
C adapted 5-Apr 2001
C
C PURPOSE:
C  1)DISTRIBUTE INACTIVE FCX AND ACTIVE FVX FOCK MATRICES AND QAX AND
C    QBX MATRICES INTO ORBITAL PART OF LINEAR TRANSFORMED VECTOR
C
C                  ( [ E(K,L) , H ] )    K<L
C       E[2]*N = - (                )
C                  ( [ E(L,K) , H ] )    K>L
C    X MAY REFER TO EITHER ONE INDEX TRANSFORMED MATRICES (N IS A
C    ORBITAL TRIAL VECTOR ) OR TRANSITION DENSITY MATRICES (N IS A
C    CONFIGURATION TRIAL VECTOR). FOR CONFIGURATION TRIAL VECTORS
C    FCX IS IDENTICALLY ZERO. FURTHER OVERLAP IS ZERO BECAUSE TRIAL
C    VECTORS ARE CHOSEN ORTHOGONAL TO REFERENCE STATE.
C
C  2)CREATE LINEAR TRANSFORMED VECTOR S[2]*N FOR N EQUAL TO EITHER A
C    ORBITAL OR A CONFIGURATION TRIAL VECTOR
C
C    ONEIND = .TRUE. FOR A ORBITAL TRIAL VECTOR
C
C                ******************************
C
C    [ E(P,Q) , H ] =
C
C    ACTIVE-INACTIVE (P,Q) = (T,I)
C                             L,K
C    FCX(I,X)*DV(T,X)-2*FCX(I,T)-2FVX(I,T)+QBX(I,T)
C        K       L          K,L       K,L      K,L
C
C    INACTIVE-ACTIVE         (I,U)
C                             K,L
C    -FCX(X,I)*DV(X,U)+2*FCX(U,I)+2*FVX(U,I)-QAX(I,U)
C           K       L        L,K        L,K      K,L
C
C    ACTIVE-ACTIVE           (T,U)
C                             K,L
C                             L,K
C    FCX(U,X)*DV(T,X)-FCX(X,T)*DV(X,U)+QBX(U,T)-QAX(T,U)
C        L       K          K       L      L,K      K,L
C        K       L          L       K      K,L      L,K
C
C    ACTIVE-SECUNDARY        (T,A)
C                             K,L
C    FCX(A,X)*DV(T,X)+QBX(A,T)
C        L       K        L,K
C
C    SECUNDARY-ACTIVE        (A,U)
C                             L,K
C    -FCX(X,A)*DV(X,U)-QAX(A,U)
C           L       K      L,K
C
C    INACTIVE-SECUNDARY      (I,A)
C                             K,L
C    2*FCX(A,I)+2FVX(A,I)
C          L,K       L,K
C
C    SECUNDARY-INACTIVE      (A,J)
C                             L,K
C    -2*FCX(J,A)-2*FVX(J,A)
C           K,L        K,L
C
C
#include "implicit.h"
#include "priunit.h"
C
      LOGICAL ONEIND, TRPDEN
C
      DIMENSION FC(*),FCX(NORBT,NORBT,*),FVX(NORBT,NORBT,*)
      DIMENSION QAX(NORBT,NASHDI,*),QBX(NORBT,NASHDI,*)
      DIMENSION UDV(NASHDI,NASHDI,*), EVECS(KZYVAR,*)
C
C  INFDIM : NASHDI
C
#include "maxorb.h"
#include "maxash.h"
#include "infvar.h"
#include "inforb.h"
#include "infind.h"
#include "infdim.h"
#include "infpri.h"
#include "infrsp.h"
#include "wrkrsp.h"
C
C -- local constants
C
      PARAMETER ( D0 = 0.0D0 , D2 =2.0D0 , DM1 = -1.0D0 )
C
C
      KYCONF = KZCONF + KZVAR
      IF (TRPDEN) THEN
         DC = D0
      ELSE
         DC = D2
      END IF
C
C DISTRIBUTE FOCK MATRICES IN EVECS
C
         KSYM1 = 0
         DO 1300 IG = 1,KZWOPT
            K     = JWOP(1,IG)
            L     = JWOP(2,IG)
            KSYM  = ISMO(K)
            LSYM  = ISMO(L)
            IF( KSYM.NE.KSYM1 ) THEN
               KSYM1 = KSYM
               NORBK = NORB(KSYM)
               IORBK = IORB(KSYM)
               NASHK = NASH(KSYM)
               NISHK = NISH(KSYM)
               IASHK = IASH(KSYM)
               IIORBK= IIORB(KSYM)
               IORBL = IORB(LSYM)
               NASHL = NASH(LSYM)
               NISHL = NISH(LSYM)
               NORBL = NORB(LSYM)
               IASHL = IASH(LSYM)
               IIORBL= IIORB(LSYM)
            END IF
            NK    = K - IORBK
            NL    = L - IORBL
            ITYPK = IOBTYP(K)
            ITYPL = IOBTYP(L)
            IF ( ITYPK.EQ.JTINAC )THEN
               DO 2000 ISIM = 1 ,NSIM
                  IF (ONEIND) THEN
                     EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                       + DC * FCX(L,K,ISIM)
                     EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                       - DC * FCX(K,L,ISIM)
                     IF ( NASHT . GT . 0 ) THEN
                        EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                       + D2 * FVX(L,K,ISIM)
                        EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                       - D2 * FVX(K,L,ISIM)
                     END IF
                  ELSE
                     EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                       + D2 * FVX(L,K,ISIM)
                     EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                       - D2 * FVX(K,L,ISIM)
                  END IF
                  IF ( ITYPL.EQ.JTACT ) THEN
                     NWL = ISW(L) - NISHT
                     IF (ONEIND) THEN
                        EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                  - DDOT(NASHL,FCX(IORBL+NISHL+1,K,ISIM),1,
     *                                        UDV(IASHL+1,NWL,1),1)
                        DO 730 IX = 1,NASHL
                           EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                     + FCX(K,IORBL+NISHL+IX,ISIM)*
     *                                        UDV(NWL,IASHL+IX,1)
 730                    CONTINUE
                     ELSE
                        TEMP1 = D0
                        TEMP2 = D0
                        NX  = NISHK
                        DO 825 NWX = IASHK+1,IASHK+NASHK
                           NX = NX + 1
                           IF (NX.LE.NK) THEN
                              FCXK = FC(IIORBK+IROW(NK)+NX)
                           ELSE
                              FCXK = FC(IIORBK+IROW(NX)+NK)
                           END IF
                           TEMP1 = TEMP1 + FCXK * UDV(NWX,NWL,ISIM)
                           TEMP2 = TEMP2 + FCXK * UDV(NWL,NWX,ISIM)
 825                    CONTINUE
                        EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                                          - TEMP1
                        EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                                          + TEMP2
                     END IF
                     EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                     - QAX(K,NWL,ISIM)
                     EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                     + QBX(K,NWL,ISIM)
                  END IF
 2000          CONTINUE
            ELSE
              IF (ITYPL.EQ.JTACT) THEN
                  NWL = ISW(L) - NISHT
                  NWK = ISW(K) - NISHT
                  DO 2020 ISIM=1,NSIM
                     IF (ONEIND) THEN
                        EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                  -DDOT(NASHL,FCX(IORBL+NISHL+1,K,ISIM),1,
     *                                        UDV(IASHL+1,NWL,1),1)
                        DO 740 IX = 1,NASHL
                           EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                        +FCX(K,IORBL+NISHL+IX,ISIM)*
     *                                        UDV(NWL,IASHL+IX,1)
 740                    CONTINUE
                        EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                     -DDOT(NASHK,FCX(IORBK+NISHK+1,L,ISIM),1,
     *                                        UDV(IASHK+1,NWK,1),1)
                        DO 750 IX = 1,NASHK
                           EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                        +FCX(L,IORBK+NISHK+IX,ISIM)*
     *                                        UDV(NWK,IASHK+IX,1)
 750                    CONTINUE
                     ELSE
                        TEMP1 = D0
                        TEMP2 = D0
                        NX  = NISHL
                        DO 835 NWX = IASHL+1,IASHL+NASHL
                           NX = NX + 1
                           IF (NX.LE.NL) THEN
                              FCXL = FC(IIORBL+IROW(NL)+NX)
                           ELSE
                              FCXL = FC(IIORBL+IROW(NX)+NL)
                           END IF
                           TEMP1 = TEMP1 + FCXL * UDV(NWX,NWK,ISIM)
                           TEMP2 = TEMP2 + FCXL * UDV(NWK,NWX,ISIM)
 835                    CONTINUE
                        EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                                          + TEMP2
                        EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                                          - TEMP1
                        TEMP1 = D0
                        TEMP2 = D0
                        NX  = NISHK
                        DO 845 NWX = IASHK+1,IASHK+NASHK
                           NX = NX + 1
                           IF (NX.LE.NK) THEN
                              FCXK = FC(IIORBK+IROW(NK)+NX)
                           ELSE
                              FCXK = FC(IIORBK+IROW(NX)+NK)
                           END IF
                           TEMP1 = TEMP1 + FCXK * UDV(NWX,NWL,ISIM)
                           TEMP2 = TEMP2 + FCXK * UDV(NWL,NWX,ISIM)
 845                    CONTINUE
                        EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                                          - TEMP1
                        EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                                          + TEMP2
                     END IF
                     EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                  +QBX(L,NWK,ISIM) -QAX(K,NWL,ISIM)
                     EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                  +QBX(K,NWL,ISIM) - QAX(L,NWK,ISIM)
 2020             CONTINUE
               ELSE
                  NWK = ISW(K) - NISHT
                  DO 2030 ISIM=1,NSIM
                     IF (ONEIND) THEN
                        EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                     -DDOT(NASHK,FCX(IORBK+NISHK+1,L,ISIM),1,
     *                                        UDV(IASHK+1,NWK,1),1)
                        DO 760 IX = 1,NASHK
                           EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                        +FCX(L,IORBK+NISHK+IX,ISIM)*
     *                                        UDV(NWK,IASHK+IX,1)
 760                    CONTINUE
                     ELSE
                        TEMP1 = D0
                        TEMP2 = D0
                        NX  = NISHL
                        DO 860 NWX = IASHL+1,IASHL+NASHL
                           NX = NX + 1
                           IF (NX.LE.NL) THEN
                              FCXL = FC(IIORBL+IROW(NL)+NX)
                           ELSE
                              FCXL = FC(IIORBL+IROW(NX)+NL)
                           END IF
                           TEMP2 = TEMP2 + FCXL * UDV(NWX,NWK,ISIM)
                           TEMP1 = TEMP1 + FCXL * UDV(NWK,NWX,ISIM)
 860                    CONTINUE
                        EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                                          + TEMP1
                        EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                                          - TEMP2
                     END IF
                     EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM)
     *                     -QAX(L,NWK,ISIM)
                     EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM)
     *                     +QBX(L,NWK,ISIM)
 2030             CONTINUE
               ENDIF
            ENDIF
 1300    CONTINUE
C
C CHANGE SIGN ON ORBITAL PART OF LINEAR TRANSFORMATION E[2]*N
C
      DO 3000 ISIM = 1,NSIM
         CALL DSCAL(KZWOPT,DM1,EVECS(KZCONF+1,ISIM),1)
         CALL DSCAL(KZWOPT,DM1,EVECS(KYCONF+1,ISIM),1)
 3000 CONTINUE
C
      IF (IPRRSP.GT.110) THEN
         WRITE(LUPRI,*)' HSOORB: LINEAR TRANSFORMATION WITH  E(2)'
         CALL OUTPUT(EVECS,1,KZYVAR,1,NSIM,KZYVAR,NSIM,1,LUPRI)
      END IF
C
C END OF HSOORB
C
      RETURN
      END
C  /* Deck hsoal2 */
      SUBROUTINE HSOAL2 (WORD,GP,CMO,UDV,PV,XINDX,MJWOP,WRK,LWRK)
C
C Alternative routine for constructing spin-orbit gradient
C by way of X2HSO Case 2
C
#include "implicit.h"
      DIMENSION GP(*), CMO(*),UDV(*),PV(*),XINDX(*),WRK(*)
      CHARACTER*8 WORD
C
      PARAMETER (D1=1.0D0, DM1=-1.0D0, DH=0.5D0)
      LOGICAL DO1,DO2,LORB,LCON,LREFST,TRPSAV
      CHARACTER*8 LABEL
C
C Used from common blocks
C
C INFINP: FLAG()
C INFORB: NORBT...
C INFVAR: NCONF
C WRKRSP: KZYVAR...,KSYMOP
C INFRSP: NCREF,IREFSY,DIROIT,SOPPA
C INFHSO: DOSO1
C INFTAP: LUINDX, LUMHSO
C
#include "priunit.h"
#include "maxorb.h"
#include "infinp.h"
#include "inforb.h"
#include "wrkrsp.h"
#include "infrsp.h"
#include "infvar.h"
#include "qrinf.h"
      DIMENSION MJWOP(2,MAXWOP,8)
#include "infhso.h"
#include "trhso.h"
#include "inftap.h"
C
      CALL QENTER('HSOAL2')
      CALL HEADER('HSOAL2: TEST OF SPIN-ORBIT GRADIENT',1)
      IF (SOPPA) CALL QUIT('HSOAL2: SOPPA not implemented yet!')
      IF (.NOT.TDHF .AND. .NOT.FLAG(27))
     &   CALL QUIT('HSOAL2 is only implemented for .DETERMINANTS')
      IF (.NOT.TDHF .AND. NCONF.NE.NCREF)
     &   CALL QUIT('HSOAL2 inconsistency : NCONF.ne.NCREF')
C
C Initialise MZYVAR... for symmetries KSYMOP and 1 (from INFVAR)
C
C Symmetry KSYMOP:
C
      CALL SETZY(MJWOP)
C
C Symmetry 1
C
      LUINDX = -1
      CALL GPOPEN(LUINDX,'LUINDF','UNKNOWN',' ','UNFORMATTED',IDUMMY,
     &            .FALSE.)
      REWIND LUINDX
      CALL MOLLAB('EXOPSYM1',LUINDX,LUERR)
      READ (LUINDX) IWOPT,IWOP,IWOPH
      CALL GPCLOSE(LUINDX,'KEEP')
      IVAR = IWOPT + NCONF
      MZVAR(1)  = IVAR
      MZCONF(1) = NCONF
      MZWOPT(1) = IWOPT
      MZYVAR(1) = 2*IVAR
      MZYCON(1) = 2*NCONF
      MZYWOP(1) = 2*IWOPT
      WRITE(LUPRI,'(/A)')    ' Number of variables symmetry 1'
      WRITE(LUPRI,'( A,I5)') '  Rotations:     ',MZWOPT(1)
      WRITE(LUPRI,'( A,I5)') '  Configurations:',MZCONF(1)
      WRITE(LUPRI,'(/A,I2)') ' Number of variables symmetry',KSYMOP
      WRITE(LUPRI,'( A,I5)') '  Rotations:     ',MZWOPT(KSYMOP)
      WRITE(LUPRI,'( A,I5)') '  Configurations:',MZCONF(KSYMOP)
      WRITE(LUPRI,'(/A)')    '  Rotational indices     K    L'
      DO 30 I=1,MZWOPT(KSYMOP)
30       WRITE(LUPRI,'(3I5)') I,MJWOP(1,I,KSYMOP),MJWOP(2,I,KSYMOP)
C
C Allocate work space
C
      KFREE = 1
      LFREE = LWRK
      CALL MEMGET('REAL',KVEC,2*IVAR,WRK,KFREE,LFREE)
      CALL MEMGET('REAL',KFI,N2ORBX,WRK,KFREE,LFREE)
C     KFI  = KVEC + 2*IVAR
C     KFREE = KFI  + N2ORBX
C     LFREE = LWRK - KFREE + 1
C     IF (LFREE.LT.0) CALL ERRWRK('HSOAL2',KFREE-1,LWRK)
C
C Convention: X1SPNORB: one-electron parts only
C             X2SPNORB: two-electron parts only
C             X SPNORB: both (default)
C                       overridden by SO1ONLY or SO2ONLY in HSOINP
C
      IF (WORD(2:2).EQ.'1' .OR. WORD(2:2).EQ.' '.AND.DOSO1) THEN
         LABEL = WORD
         LABEL(2:2) = '1'
         KSYMP = -1
         CALL PRPGET(LABEL,CMO,WRK(KFI),KSYMP,ANTSYM,
     &               WRK(KFREE),LFREE,IPRHSO)
         IF (KSYMP.NE.KSYMOP) THEN
            WRITE (LUPRI,'(/A/2A/A,2I5/A,F10.2)')
     &           'FATAL ERROR: KSYMOP .ne. KSYMP from PRPGET',
     &           '   Property label  : ',LABEL,
     &           '   KSYMOP and KSYMP:',KSYMOP,KSYMP,
     &           '   ANTSYM          :',ANTSYM
            CALL QUIT('KSYMOP .ne. KSYMP from PRPGET')
         END IF
      ELSE
         CALL DZERO(WRK(KFI),N2ORBX)
      END IF
C
C
C Set up the vector to be  N_o = ( 0 ) , N_c = ( - Ref )
C                                ( 0 )         (   Ref )
C
      CALL DZERO(WRK(KVEC),2*IVAR)
      CALL GETREF(WRK(KVEC + IVAR),NCONF)
      CALL DAXPY(NCONF,DM1,WRK(KVEC + IVAR),1,WRK(KVEC),1)
      WRITE(LUPRI,'(/A)') 'REFERENCE CI  VECTOR'
      CALL OUTPUT(WRK(KVEC + IVAR),1,NCONF,1,1,NCONF,1,1,LUPRI)
      WRITE(LUPRI,'(/A)') 'CALLING HSO2CR WITH VECTOR'
      CALL RSPPRC(WRK(KVEC),NCONF,IVAR,LUPRI)
      CALL RSPPRW(WRK(KVEC+NCONF),MJWOP,IWOPT,1,IVAR,LUPRI)
      ISYMV = 1
      ISPINV = 0
      IKLVL = 0
      DIROIT = .TRUE.
C
C Gradient setup
C
      IGRSYM = KSYMOP
      IF (TRPLET) THEN
         IGRSPI = 1
      ELSE
         IGRSPI = 0
      END IF
C
C Operator setup
C
      IOPSYM = KSYMOP
      IOPSPI = 1
C
C Orbital gradient
C
      LORB = KZWOPT.GT.0
      LPV = 2*N2ASHX*N2ASHX
      CALL MEMGET('REAL',KD,N2ASHX,WRK,KFREE,LFREE)
      CALL MEMGET('REAL',KP,LPV,WRK,KFREE,LFREE)
      KP1=KP
      KP2=KP+N2ASHX*N2ASHX
      IF (TRPLET) THEN
         CALL DCOPY(N2ASHX,UDV,1,WRK(KD),1)
         CALL DCOPY(LPV,PV,1,WRK(KP),1)
         CALL IPSET(0,0,1,1)
      ELSE
         IF (TDHF) THEN
            IF (NASHT.EQ.1) THEN
               WRK(KD)=1.0D0
            END IF
         ELSE
            CALL MEMGET('REAL',KCREF,MZCONF(1),WRK,KFREE,LFREE)
            CALL GETREF(WRK(KCREF),MZCONF(1))
            CALL RSPDM(IREFSY,IREFSY,MZCONF(1),MZCONF(1),
     &      WRK(KCREF),WRK(KCREF),
     &      WRK(KD),WRK(KP),0,1,.FALSE.,.FALSE.,XINDX,
     &      WRK,KFREE,LFREE)
            CALL MEMREL('HSOCTL<-RSPDM',WRK,KD,KCREF,KFREE,LFREE)
            CALL MTRSP(N2ASHX,N2ASHX,WRK(KP1),N2ASHX,WRK(KP2),N2ASHX)
         END IF
         CALL IPSET(0,1,1,0)
      END IF
      OVLAP = D1
      ISYMDN = 1
C
C Configurational gradient
C
      IF ( IGRSYM.EQ.1) THEN
         LCON = KZCONF.GT.1
      ELSE
         LCON = KZCONF.GT.0
      END IF
      LREFST = .FALSE.
C
C Create the gradient
C
      TRPSAV=TRPLET
      CALL HSO2CR(IGRSYM,IGRSPI,GP,WRK(KVEC),
     *            2*IVAR,NCONF,ISYMV,ISYMDN,
     *            XINDX,OVLAP,WRK(KD),WRK(KP),WRK(KFI),
     *            WRK(KFREE),LFREE,
     *            KZYVAR,LCON,LORB,LREFST,LUMHSO,KSYMOP,
     *            IOPSPI,ISPINV,IKLVL,DUM,IDUM,DUM,IDUM,MJWOP)
      CALL HEADER('SPIN-ORBIT GRADIENT FROM HSOAL2',3)
      TRPLET=TRPSAV
C
      IF (IPRRSP.GT.100) THEN
         WRITE(LUPRI,'(/A)') ' Orbital part'
         CALL OUTPUT(GP(1+KZCONF),1,KZWOPT,1,2,KZVAR,2,1,LUPRI)
         WRITE(LUPRI,'(/A)') ' Configuration part'
         CALL OUTPUT(GP,1,KZCONF,1,2,KZVAR,2,1,LUPRI)
      ELSE
         CALL RSPPRW(GP(1+KZCONF),MJWOP,KZWOPT,KSYMOP,KZVAR,LUPRI)
         CALL RSPPRC(GP,KZCONF,KZVAR,LUPRI)
      END IF
      IF (KSYMOP.EQ.1 .AND. LCON) THEN
         WRITE(LUPRI,'(/A)') 'Project out reference component'
         KCREF = KVEC + IVAR
         OVLAP = DDOT(KZCONF,WRK(KCREF),1,GP,1)
         CALL DAXPY(KZCONF,-OVLAP,WRK(KCREF),1,GP,1)
         OVLAP = DDOT(KZCONF,WRK(KCREF),1,GP(1+KZVAR),1)
         CALL DAXPY(KZCONF,-OVLAP,WRK(KCREF),1,GP(1+KZVAR),1)
         CALL RSPPRW(GP(1+KZCONF),MJWOP,KZWOPT,KSYMOP,KZVAR,LUPRI)
         CALL RSPPRC(GP,KZCONF,KZVAR,LUPRI)
      END IF
      CALL MEMREL('HSOAL2',WRK,1,1,KFREE,LFREE)
      CALL QEXIT('HSOAL2')
      RETURN
      END
C
C /* Deck get_FSO_AO */
C
      SUBROUTINE GET_FSO_AO(WORD,TRIPLET,F,D,ND)
C
C Generalize GETFMAT to handle non-symmetric densities (quadratic
C response spin-orbit)
C
#include "implicit.h"
      CHARACTER*8 WORD
      INTEGER ND
      LOGICAL TRIPLET(ND)
#include "inforb.h"
      DIMENSION F(NBAST,NBAST,ND), D(NBAST,NBAST,ND)

#include "infhso.h"
#include "priunit.h"
#include "iratdef.h"
#include "eribuf.h"
C
C LOCAL
C
      PARAMETER (LBUF_alloc = 600)
      REAL*8    BUF(2*LBUF_alloc+1)
      INTEGER   IINDX4(4,LBUF_alloc), LENGTH
      INTEGER   IB, I,J,K,L,ICOOR, LUHSOAO, COMPONENT,ID
      REAL*8    LEFT(2), RIGHT(2)
      REAL*8    GIJKL,XIJKL,DIJ,DKL, DP5, D1P5, D1, D2 
      PARAMETER(DP5=0.5D0, D1P5=1.5D0, D1=1.0D0, D2=2.0D0)
C
      CALL QENTER('GET_FSO_AO')
      IF (ND .GT. 2) THEN
         CALL QUIT('ND .gt. fixed dimension of 2')
      END IF
C
C READ AND DISTRIBUTE SPIN-ORBIT AO INTEGRALS
C
      DO ID=1,ND
         IF (TRIPLET(ID)) THEN
            LEFT(ID)= D1
            RIGHT(ID)=D2
         ELSE
            LEFT(ID) =D2
            RIGHT(ID)=D1
         END IF
      END DO

      LUHSOAO=-1
      CALL GPOPEN(
     &   LUHSOAO,'AO2SOINT','OLD','SEQUENTIAL','UNFORMATTED',-1,.FALSE.
     &   )
      REWIND (LUHSOAO)

      CALL ERIBUF_INI  ! set NIBUF, NBITS, IBIT1, IBIT2
      LBUF = LBUF_alloc

      LEN_BUF4 = 2*LBUF + NIBUF*LBUF + 1 ! BUF(LBUF), IBUF4(NIBUF,LBUF), LENGTH4 in integer*4 units

      CALL MOLLAB('AO2SOINT',LUHSOAO,LUPRI)
      IF (WORD(1:1).EQ.'X') COMPONENT = 1
      IF (WORD(1:1).EQ.'Y') COMPONENT = 2
      IF (WORD(1:1).EQ.'Z') COMPONENT = 3
      
      IF (IPRHSO .GT. 100 .OR. HSODEBUG) THEN
         WRITE(LUPRI,*) 'Hello from GET_FSO_AO. WORD: ',WORD,COMPONENT
         WRITE(LUPRI,*) ' TRIPLET? ',TRIPLET(1:ND)
         WRITE(LUPRI,*) ' LBUF, NIBUF, NBITS, IBIT1, IBIT2, LEN_BUF4',
     &                    LBUF, NIBUF, NBITS, IBIT1, IBIT2, LEN_BUF4
      END IF

      NRECS = 0
cF90  READFILE: DO 
 100  CONTINUE
         CALL READI4(LUHSOAO, LEN_BUF4, BUF)
         NRECS = NRECS + 1
         CALL AOLAB4(BUF(LBUF+1),LBUF,NIBUF,NBITS,IINDX4,LENGTH)
         IF (HSODEBUG) THEN
         IF (NRECS .LE. 2) THEN
            write (lupri,*) 'AO2SOINT record no., length:',NRECS, length
            do ib = 1,length
               write(lupri,*) iindx4(1:4,ib), buf(ib)
            end do
         END IF
         END IF
cF90     IF (LENGTH.LE.0) EXIT
         IF (LENGTH.LE.0) GO TO 300
cF90     BUFFER: DO IB=1,LENGTH
         DO IB=1,LENGTH
            K = IINDX4(3,IB)
            L = IINDX4(4,IB)
            IF (L.EQ.0) THEN
               ICOOR = K
cF90           CYCLE BUFFER
               GO TO 200
            END IF
cF90        IF (ICOOR.NE.COMPONENT) CYCLE BUFFER
            IF (ICOOR.NE.COMPONENT) GO TO 200
            I = IINDX4(1,IB)
            J = IINDX4(2,IB)
            GIJKL=BUF(IB)
            IF (I.EQ.J) GIJKL = DP5*GIJKL
            XIJKL=D1P5*GIJKL
            DO ID=1,ND
               DIJ = LEFT(ID)*(D(I,J,ID) + D(J,I,ID))
               DKL = RIGHT(ID)*(D(L,K,ID) - D(K,L,ID))
               F(I,J,ID) = F(I,J,ID) + GIJKL*DKL
               F(J,I,ID) = F(J,I,ID) + GIJKL*DKL
               F(K,L,ID) = F(K,L,ID) + GIJKL*DIJ
               F(L,K,ID) = F(L,K,ID) - GIJKL*DIJ
               F(I,L,ID) = F(I,L,ID) - XIJKL*D(J,K,ID)
               F(L,I,ID) = F(L,I,ID) + XIJKL*D(K,J,ID)
               F(I,K,ID) = F(I,K,ID) + XIJKL*D(J,L,ID)
               F(K,I,ID) = F(K,I,ID) - XIJKL*D(L,J,ID)
               F(J,L,ID) = F(J,L,ID) - XIJKL*D(I,K,ID)
               F(L,J,ID) = F(L,J,ID) + XIJKL*D(K,I,ID)
               F(J,K,ID) = F(J,K,ID) + XIJKL*D(I,L,ID)
               F(K,J,ID) = F(K,J,ID) - XIJKL*D(L,I,ID)
            END DO
cF90     END DO BUFFER
  200    CONTINUE
         END DO
cF90  END DO READFILE
      GO TO 100
  300 CONTINUE
C
      CALL GPCLOSE(LUHSOAO,'KEEP')
      CALL QEXIT('GET_FSO_AO')
cF90  END SUBROUTINE GET_FSO_AO
      END
